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Preface 


This  study  is  another  link  in  a  growing  chain  of  research  at  the 
Air  Force  Institute  of  Technology  to  design  a  tracking  algorithm  for 
use  with  the  Air  Force  Weapons  Laboratory's  high  energy  laser  weapon 
system.  As  such,  my  work  extends  the  multiple  model  adaptive  filter 
developed  by  Lt.  Robert  I.  Suizu  and  investigates  Bayesian  vs  maximum 
a  posteriori  estimation  using  the  multiple  model  adaptive  filter 
concept. 

I  owe  alot  to  my  predecessors  in  this  effort;  without  whom  an 
investigation  of  this  complexity  could  not  have  been  attempted.  My 
sincere  gratitude  goes  to  Dr.  Peter  S.  Maybeck,  my  thesis  advisor, 
for  the  effort  and  time  expended  on  my  behalf.  His  ability  to 
consistently  fan  the  flames  of  my  enthusiasm  and  challenge  my  mind 
contributed  significantly  to  keeping  my  work  on  track  and  productive. 
Finally,  I  thank  my  friends,  classmates,  and  instructors  whose  sense 
of  humor  was  there  when  mine  was  not. 
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Abstract 


Previous  efforts  at  the  Air  Force  Institute  of  Technology  have 
led  to  the  development  of  a  multiple  model  adaptive  filter  (MMAF) 
tracking  algorithm  which  demonstrates  significant  improvements  in 
performance  against  close-range,  highly  dynamic,  airborne  targets, 
over  a  direct  correlation  method  currently  in  use.  The  basic 
elemental  filter  in  the  MMAF  bank  combines  an  enhanced  correlator  and 
a  linear  Kalman  filter.  Digital  signal  processing  techniques  are 
used  to  derive  a  target  shape  function  from  the  forward  looking 
infrared  sensor  data.  This  shape  function  is  used  as  a  template  in 
the  correlation  algorithm  which  generates  offset  pseudo-measurements 
for  the  update  portion  of  a  linear  Kalman  filter.  The  multiple 
models  are  created  by  tuning  the  basic  model  for  "best"  performance 
against  differing  target  maneuvering  behavior  and  with  physically 
different  fields  of  view.  The  outputs  of  three  independent  elemental 
filters,  each  receiving  data  from  a  shared  sensor  are  used  to 
generate  a  single  adaptive  estimate  of  the  state  via  a  probabilistic 
weighted  average  (Bayesian  form)  or  by  selection  of  the  one  elemental 
filter  associated  with  the  highest  probability  (MAP  form).  The 
adaptive  state  estimate  can  produce  target  position  predictions  to  be 
used  in  generating  feedback  control  for  maintaining  the  target  in  the 
center  of  the  field  of  view. 

There  are  two  main  results  from  this  effort.  The  addition  of  a 
third  elemental  filter  to  the  baseline  MMAF  improves  tracking  per¬ 
formance  over  the  two-element  MMAF.  Specifically,  the  peak  error 
following  a  maneuver  is  significantly  reduced.  However,  the  MAP 
estimation  approach  does  not  differ  significantly  from  the  Bayesian 


approach . 


BAYESIAN  VS  MAP 

MULTIPLE  MODEL  ADAPTIVE  ESTIMATION 
FOR  FIELD  OF  VIEW  EXPANSION  IN 
TRACKING  AIRBORNE  TARGETS 


I.  INTRODUCTION 

Since  the  inception  of  the  Strategic  Defense  Initiative  (SDI) , 
increasing  interest  has  been  focused  on  laser  weapon  systems.  The 
extremely  tight  specif ications  required  for  the  pointing  and  tracking 
systems  supporting  these  laser  weapons  has  motivated  research  into 
innovative  methods  of  accurately  tracking  targets  maneuvering  at  high 
velocities . 

1 . 1  Background 

The  Air  Force  Weapons  Laboratory  at  Kirtland  AFB,  New  Mexico, 
maintains  ongoing  research  efforts  using  high  energy  laser  weapons 
against  airborne  targets.  Currently,  pointing  and  tracking  functions 
are  accomplished  by  a  correlation  algorithm.  A  Forward  Looking  Infra- 
Red  (FLIR)  sensor  makes  target  measurements,  which  the  operational 
correlator  uses,  comparing  data  from  the  previous  sample  period  to  data 
available  during  the  current  sample  time.  Cross  correlations  of  the 
data  sets  generate  relative  positions  offsets,  presumably  due  to  target 
motion.  These  offsets  are  used  to  center  the  target  in  the  FLIR  field 
of  view  (FOV) .  Allowing  for  boresight  error  corrections,  the  laser  is 
coupled  to  the  FLIR  so  that  positioning  the  center  of  the  sensor  FOV 
also  points  the  laser  toward  the  target. 


Although  the  correlation  tracker  performs  reasonably  well  against  a 
wide  variety  of  targets,  it  has  several  inherent  limitations.  The 
algorithm  processes  the  most  recent  data,  providing  information  on  the 
current  position  of  the  target  relative  to  the  target  position  at  the 
last  sample  time.  Controls  are  applied  some  finite  calculation  time 
later.  This  delay  causes  pointing  errors,  which  become  increasingly 
detrimental  with  highly  maneuverable  targets.  Additionally,  the 
relative  offsets  are  a  function  of  image  changes  between  one  frame  and 
the  next,  deviations  not  entirely  due  to  target  motion.  These 
deviations  can  also  be  attributed  to  mirror  vibration,  atmospheric 
distortion,  and  inherent  FLIR  measurement  errors.  The  correlator  does 
not  differentiate  between  the  effects  of  these  error  sources  and  actual 
target  motion.  Research  into  alternate  tracking  algorithms  is  motivated 
by  the  above  limitations .  Since  1978,  the  Air  Force  Institute  of 
Technology  has  supported  a  number  of  Master's  theses  demonstrating  the 
feasibility  and  performance  of  a  tracking  algorithm  based  on  Kalman 
filtering  techniques. 

The  Kalman  filter  directly  addresses  the  previously  mentioned 
limitations,  by  attempting  to  separate  actual  target  motion  from  other 
phenomena  that  are  observed  as  apparent  target  motion.  Using  the 
statistical  characteristics  of  atmospheric  jitter  and  measurement 
errors,  and  a  model  of  the  anticipated  target  dynamics,  the  filter 
predicts  an  estimated  target  position.  This  prediction  allows  the 
FLIR  to  anticipate  target  motion,  thus  reducing  tracking  error  due  to 


time  delays. 


The  initial  feasibility  study  by  Mercier  [6],  based  on  a  four- 
state  extended  Kalman  filter,  assumed  the  distant  target's  intensity 
pattern  on  the  FLIR  image  plane  could  be  well  modeled  as  a  bivariate 
Gaussian  distribution.  This  filter  was  based  on  a  benign  target 
dynamics  model  with  target  position  represented  by  a  zero  mean,  first 
order  Gauss-Markov  process.  The  filter's  model  for  noise  in  the 
individual  FLIR  pixel  intensity  measurements  was  uncorrelated  in  time 
and  space.  This  algorithm  produced  an  order  of  magnitude  improvement 
in  tracking  error  over  the  correlation  algorithm,  in  benign  scenarios 
with  filter  parameters  well-matched  to  the  real  world  environment. 

In  robustness  studies  conducted  by  Hamly  and  Jensen  [3],  the 
spread,  size,  and  shape  of  the  target  intensity  pattern  were  varied  to 
study  the  filter's  sensitivities  to  incorrect  assumptions.  More 
dynamic  target  motion  and  adaptive  estimation  of  target  shape,  still 
assuming  bivariate  Gaussian  profiles,  but  with  generally  noncircular 
constant  intensity  contours,  were  incorporated  and  the  filter's  dynamic 
model  was  expanded  to  estimate  velocity  and  acceleration.  Spatial  and 
temporal  noise  correlations  were  investigated,  leading  to  the  addition 
of  a  nearest-two-pixels  spatial  correlation  for  the  filter's  noise 
model.  In  an  effort  to  track  the  maneuvering  target,  gain  changes  and 
other  adaptations  based  on  ad  hoc  responses  to  a  maneuver  detection  were 
added  to  the  filter.  While  the  results  were  encouraging,  the  filter 
still  experienced  difficulty  with  harsh  maneuvers  differing  from  the 


assumed  target  dynamics  model. 


Up  until  this  point,  the  filters  have  all  assumed  a  known  target 
intensity  shape  function  with  a  single  peak.  In  contrast,  Singletery 
[12]  and  Rogers  [11]  developed  signal  processing  techniques  to  derive 
the  target  shape  function  based  on  the  intensity  measurements.  This 
technique  was  tested  against  multiple-hot-spot  targets  with  dynamic 
variations.  These  efforts,  concentrating  on  the  shape  function 
derivation,  only  pursued  the  benign  target  dynamics  (as  seen  in 
Mercier's  formulation).  Rogers  developed  an  alternative  filter 
algorithm  using  the  target  shape  function  as  a  template  for  an 
enchanced  correlator,  rather  than  using  it  as  a  measurement  function 
for  an  extended  Kalman  filter.  The  correlator  produced  offsets,  x  and 
y,  from  the  center  of  the  sensor  FOV  which  were  used  as  "measurement" 
inputs  to  a  linear  Kalman  filter.  The  linear  Kalman  filter  requires  a 
much  lower  level  of  computational  resources  than  the  extended  Kalman 
filter  and  would  be  preferable  if  comparable  performance  were  achieved. 
In  fact,  comparable  RMS  tracking  errors  were  produced  in  performance 
analyses,  with  the  extended  Kalman  filter  yielding  a  larger  mean  errors, 
while  the  linear  filter/correlator  combination  produced  larger  standard 
deviations . 

Follow-on  research  by  Millner  [8]  and  Kozemchak  [4],  incorporated 
more  accurate  target  dynamics  models.  Both  the  extended  Kalman  filter 
and  the  linear  filter/correlator  algorithm  were  tested  against  close 
range  targets,  with  maneuvers  of  up  to  20  g's.  Again,  the  extended 
Kalman  filter  exhibited  larger  biases  and  smaller  standard  deviations 
than  the  linear  filter/correlator.  However,  both  filter  formulations 


were  slow  to  respond  to  harsh  maneuvers  that  differed  significantly  from 
the  filter's  assumed  target  dynamics  models  and  the  tradeoff  between 
precision  tracking  of  benign  targets  versus  maintaining  lock  on  a  highly 
maneuvering  target  was  not  completely  resolved. 

In  an  effort  to  address  this  limitation,  Flynn  [2]  initially 
investigated  a  multiple  model  adaptive  filter  (MMAF) .  The 
attractiveness  of  this  approach  was  increased  by  the  potential  for 
distributed  processing.  The  MMAF  was  successfully  implemented  by  Suizu 
[13].  Based  on  probabilistic  weighting,  the  filter  adaptively  changed 
target  dynamics  model  and  FOV  size,  allowing  the  filter  to  maintain  lock 
and  track  targets  performing  20  g  pull-ups  at  ranges  of  20  km.  This 
approach  was  tested  using  both  the  extended  Kalman  filter  and  the  linear 
filter/correlator  with  both  exponentially  time-correlated  acceleration 
and  constant  turn-rate  dynamics  models.  The  two  types  of  filters 
continued  the  same  relative  performance  as  seen  in  [8]  and  [4],  but  were 
capable  of  tracking  targets  with  trajectories  that  differed 
significantly  in  maneuver  magnitude  during  the  course  of  the  tracking 
scenario. 

1 .2  Problem 

This  effort  concentrates  on  expanding  the  multiple  model 
linear/correlator  algorithm  developed  by  Suizu  [13].  The  potential  for 
decreased  computational  loading’  compared  to  the  extended  Kalman  filter, 
while  maintaining  comparable  accuracy,  makes  this  filter  more  attractive 
for  further  development.  The  model,  as  developed,  will  be  operated  at  a 
100  Hz  sampling  rate,  as  well  as  30  Hz,  to  investigate  the  potential 
benefit  of  optical  implementation  of  the  algorithms  discussed  by  Rogers 


[11]  and  Roemer  [10].  Additionally,  the  simulation  will  be  extended 
beyond  the  current  five  seconds;  both  filters  currently  exhibit  some 
difficulty  with  the  minimum  range/maximum  crossing  rate  conditions  at 
the  5  second  point  in  the  simulation,  and  this  will  allow  exploration  of 
tracking  performance  beyond  this  difficult  period.  To  improve  overall 
performance  against  highly  maneuvering  targets,  another  filter  is 
incorporated  in  the  MMAF  tracker.  Additionally,  Maximum  a  Posteriori 
(MAP)  estimation  is  compared  to  Bayesian  estimation,  i.e.,  producing  an 
adaptive  state  estimate  from  the  one  model  associated  with  the  highest 
probability  of  validity  rather  than  from  a  probabilistically  weighted 
average  of  all  models.  This  comparison  explores  the  tradeoff  between 
faster  response  to  maneuvers  and  decreased  accuracy  of  the  estimate. 
Within  the  MMAF  is  a  bank  of  N  separate  filters.  Each  of  these  can  be 
an  extended  Kalman  filter  or  a  linear  filter /correlator  combination  in 
this  problem. 

1.2.1  The  Extended  Kalman  Filter  Tracker.  The  extended  Kalman 
filter  algorithm  uses  measurements  from  the  FL1R  sensor  to  update  the 
state  estimate,  which  is  propagated  forward  in  time,  based  on  the 
dynamics  model  of  the  Kalman  filter  and  used  to  control  the  FLIR/ laser 
orientation  as  shown  in  Figure  1-1.  The  measurement  vector  £(t^) 
consists  of  the  intensity  output  of  64  pixels  arranged  in  an  8-by-8 
"tracking  window".  The  extended  Kalman  filter  uses  the  linearized  and 
nonlinear  intensity  functions  (  HjxCt^),^]  and  h^[jc(t^) , t^]  ),  along 

with  the  measurement  vector  £(t^),  to  update  the  state  estimate  based  on 
the  equation 


(1-1) 


xCtj^  )  =  x^)  +  Kui)[z.ui)  -  h(x(t1),ti>] 


where 

A  + 

jc(ti)  =  state  estimate  after  update 
xCt^  »  state  estimate  prior  to  update 
K(t  )  -  Kalman  filter  gain 

_z(t  )  =  measurement  vector;  the  assumed  model  is: 

£(ti)  *h[x(t1),ti]  +  v(ti) 

h^tacCt^)  ,t^]  =  intensity  shape  function  for  measurements  at  as  a 
function  of  the  state  estimate 

HfxCt^.t^]  ■  3h/3x;  the  linearized  intensity  function  evaluated  at 
jc(t-)  at  time  t-  used  to  generate  K(t^) 

v(t^)  «  measurement  noise 

Based  on  the  internal  dynamics  models,  the  Kalman  filter  will  propagate 
the  estimate  forward  one  sample  period.  The  controller  uses  this 
information  to  orient  the  FLIR  so  that  the  center  of  the  FOV  is  pointing 
at  the  predicted  target  position.  It  is  assumed  that  the  controller  is 
capable  of  repositioning  in  less  than  one  sample  period. 

The  upper  path  of  Figure  1-1  details  the  generation  of  the 
intensity  shape  functions  at  each  sample  point.  First,  the  input  data 
from  the  FLIR  is  Fourier  transformed;  the  transformation  is  motivated  by 
the  possibility  of  optical  processing  and  the  comparative  ease  of 
performing  the  necessary  computations  in  the  frequency  domain.  The 
original  8-by-8  tracking  window  is  expanded  into  a  24-by-24  data  array, 
since  the  larger  array  reduces  errors  due  to  edge  effects,  aliasing,  and 


leakage  conditions  encountered  in  transforming  a  finite  sequence 
[12:18],  In  this  application,  the  larger  24-by-24  array  is  produced  by 
padding  the  8-by-8  tracking  data  with  additional  intensity  data 
available  from  the  FLIR. 

This  transformed  24-by-24  data  array  is  altered  so  as  to  center  the 
target  intensity  pattern  and  then  temporally  averaged  with  the  previous 
frames  of  transformed  and  centered  data,  to  generate  the  desired 
nonlinear  intensity  function.  This  insures  that  any  variations  in  the 
data  are  primarily  due  to  noise  and  not  changes  in  the  target  position 
in  the  FOV.  The  "centering"  in  the  original  domain  is  actually 
performed  by  multiplying  the  transformed  data  by  a  negating  phase  shift, 
based  on  the  shifting  property  of  the  Fourier  transform  [9],  The 
negating  phase  shift  is  the  complex  conjugate  of  the  linear  phase  shift 
due  to  the  (estimated)  target  image  offset  in  the  spatial  domain.  This 
offset  is  obtained  from  the  updated  state  estimates  of  the  extended 
Kalman  filter. 

After  centering,  the  data  is  filtered  by  exponential  smoothing. 
This  exponential  smoothing  approximates  a  true  time  averaging  without 
the  corresponding  need  to  maintain  numerous  frames  of  data  in  memory. 
Because  the  noise  is  expected  to  vary  more  from  frame  to  frame  than  the 
target  intensity  pattern,  its  effect  is  reduced  by  averaging  successive 
frames  of  data.  The  resulting  image  is  used  as  the  target  intensity 
shape  function,  h[x(t^),tj.  The  linearized  intensity  function, 
H[3c(t^) ,t^] ,  is  obtained  by  applying  the  derivative  property  of  the 
Fourier  transform,  which  requires  only  a  simple  multiplication. 


These  linear  and  nonlinear  shape  functions  are  used  to  update  the 
state  estimates  at  the  next  measurement.  To  insure  the  proper  location 
of  the  shape  function  in  the  FOV,  both  functions  are  evaluated  using  the 
propagated  state  estimates,  x(t^+^).  It  is  assumed  that  the  FLIR  is 
centered  on  the  position  predicted  by  the  target  dynamics  model.  Thus, 
the  intensity  functions  are  evaluated  at  the  position  predicted  by  the 
atmospheric  states,  if  there  is  a  single  elemental  filter;  this  must  be 
modified  in  the  multiple  model  adaptive  filter,  and  the  modification  is 
discussed  in  Section  4.4.  The  Fourier  transform  shifting  property  is 
applied  again  to  phase  shift  the  transformed  image  by  the  estimated 
atmospheric  offset.  Finally,  the  inverse  Fourier  transform  is 
performed,  returning  the  intensity  functions  to  the  spatial  domain  for 
use  in  the  next  update  cycle. 

1.2.2  Linear  Kalman  Filter /Correlation  Tracker.  The  linear 
filter/correlator  developed  by  Rogers  [11]  is  very  similar  to  the 
extended  Kalman  filter  algorithm  depicted  in  Figure  1-1.  The  evaluation 
of  H  is  no  longer  required;  but  otherwise,  the  upper  path  remains  the 
same.  The  h  and  jz  are  not  used  directly  by  the  Kalman  filter.  The 
intensity  function,  h^f ac ( t ,t:i-«-l^  *  *s  usec*  as  a  template  by  an 
enhanced  correlator  algorithm.  Sensor  measurements  are  compared  to  the 
template,  producing  estimated  offsets  between  the  target  centroid  and 
the  center  of  the  sensor  FOV.  These  two  scalar  offsets  are  pseudo¬ 
measurement  inputs  to  a  linear  Kalman  filter  (the  filter  of  Section 
1.2.1  was  nonlinear  only  in  the  measurement  update  portion  of  the 
algorithm  when  a  linear  state  propagation  model  is  used) .  In  all  other 
respects,  the  two  tracking  algorithms  are  the  same. 


1 . 3  Overview 


The  following  chapters  expand  the  details  of  the  algorithm,  the 
testing  environment,  and  the  results.  Chapter  II  presents  the 
mathematical  foundation  for  the  multiple  model  algorithm  for  Bayesian 
estimation  and  MAP  estimation.  The  truth  model  is  covered  in  Chapter 
III,  establishing  the  standard  for  performance  evaluation.  Chapter  IV 
describes  the  basic  linear  filter /correlator  model  in  detail  along  with 
the  MMAF  implementation.  Chapter  V  covers  the  tracker  statistics, 
parameters,  plot  formats,  and  the  test  scenarios.  Finally,  Chapter  VI 
contains  the  performance  analysis.  The  conclusions  and  recommendations 


II.  THEORETICAL  DEVELOPMENT  OF  MULTIPLE  MODEL  ADAPTIVE  ESTIMATION 


From  a  Bayesian  point  of  view,  it  would  be  desirable  to  be  able  to 

generate  the  conditional  density  function  of  the  state  given  the 

measurement  history;  with  this  density  function,  a  good  estimate  of  the 

state  vector,  x,  is  obtained  by  evaluating  E{3^(tQ)l  Z^} ,  the  conditional 

expectation  of  the  state  given  the  measurement  history.  If  the  target 

motion  is  adequately  described  by  the  linear  Kalman  filter  dynamics 

model,  this  conditional  mean  is  easily  calculated.  However,  in  this 

problem  the  appropriate  target  dynamics  model  may  be  uncertain  and/or 

time  varying.  Let  a  be  the  vector  of  uncertain  parameters  in  the 

Kalman  filter  dynamics  model.  In  order  to  estimate  the  state  vector 

properly,  it  is  also  necessary  to  estimate  a,  thus  motivating  the 

evaluation  of  f  .  ,  the  joint  density  of  x  and  a  given  the  measurement 
x,a|4  —  — 

history.  The  evaluation  of  this  density  function  is  the  core  of 

Bayesian  estimation  [5:129]. 

2. 1  Mathematical  Development 

Application  of  Bayes  theorem  to  the  joint  conditional  density 
function  of  the  variables  to  be  estimated,  given  the  measurements, 
yields : 


fx,a|Z  =  fx|a,Z  ’  fa|Z  (2_1) 

The  first  term  on  the  right  hand  side  of  the  equal  sign  is  a  Gaussian 
distribution,  if  the  adequate  system  models  are  in  the  form  of  a  linear 
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state  differential  and  linear  measurement  equation,  all  driven  by  white 
Gaussian  noises.  As  such,  it  is  totally  specified  by  its  mean  and 


covariance,  which  are  the  outputs  of  a  Kalman  filter  based  on  one 
particular  value  of  a.  The  second  term  can  be  expressed  as 


where  a  is  a  dummy  variable  associated  with  a_  values,  and 


and 


(2-2) 


=  the  composite  vector  of  all  measurements  prior 
to  time  t. 

i 

z.  =  the  vector  of  measurements  available  at 
1  time  t^ 

The  first  numerator  term  in  Equation  (2-2)  is  a  Gaussian  density 

~  _  -  T 

function  with  mean  H(t^)x(t^)  and  covariance  [ll(t^)P^(t:^)H^  (t^)+jt(t^)  ] , 

calculated  by  a  Kalman  filter  based  on  a  given  value  of  Starting 

with  f  |  =  f  (a),  the  a  priori  probability  of  a  particular  dynamics 

a|  L  a 

model  correctly  modeling  the  actual  target  motion  at  the  initial  time 
(letting  determine  the  applicable  dynamics  model) ,  Equation  (2-2)  can 
be  calculated  recursively.  Having  evaluated  both  terms  in  Equation 
(2-1),  the  desired  conditional  mean  E(x(t^)}  ,  is  readily  determined. 


=  JS fx|z  dC 

— ao 

OO  _  OO  — t 

-  fECfx,a|2d“3  dE 

— ccr  — oo 

=  /"e{x|  a,Z  }  f  i  da 

_c/  1  al z 

This  result  is  based  on  Equation  (2-1),  the  concept  of  marginal 
densities,  Bayes'  theorem,  and  changing  the  order  of  integration.  While 
the  above  calculations  achieve  the  desired  estimate,  they  are 
computationally  unattractive.  On-line  implementation  of  the  integrals 
is  difficult  and  time  consuming,  if  at  all  tractable.  The  additional 
requirement  for  an  infinite  number  of  Kalman  filters  to  span  the 
continuous  a^  parameter  space  suggests  that  a  judicious  discretization 
might  provide  a  more  feasible  approach. 

Discretization  of  the  parameter  space  allows  the  integrals  over 
to  be  evaluated  as  the  summation  of  all  possible  a  values.  Assume  the 
parameter  £  can  take  on  a  finite  number  of  values,  selected  to  span  t'^e 
set  of  possible  values  in  the  continuous  parameter  space  and  represented 
by  the  set  . . .a^} .  Each  element  a^  represents  a  different  system 

model,  which  is  ultimately  reflected,  in  this  application,  as  a  separate 
Kalman  filter  based  on  the  dynamics  model  (and  appropriate  FOV  size) 
defined  by  a^.  Bayesian  estimation,  using  the  outputs  of  the  N  filters, 
requires  an  a  priori  density  function  for  a  as  well  as  ]t(to) .  Let 
Pk(tQ)  represent  the  probability  that  the  model  specified  by  a^  is  the 
best  representation  of  the  true  target  dynamics  at  the  initial  time. 
Also  define  the  hypothesis  conditional  probability  p.  (t.)  as 


The  recursion  for  p^t^)  is  expressed  in  terms  of  p^t^  ^)  and  other 
elements  that  can  be  evaluated.  In  terms  of  these  p^t^)  values. 
Equation  (2-3)  becomes 


x(t£)  = 


k=l  pk(ti) 

k-1  *k(tt>  Pk(ti) 


(2-6) 


where  x^(t*)  is  produced  by  the  Kalman  filter  that  assumes  the 
parameter  vector  equals  a^.  The  denominator  of  Equation  (2-5)  is  a 
normalizing  constant  that  is  the  same  for  each  conditional  probability, 
p  (t.),  and  is  the  sum  of  all  the  numerator  terms.  If  these 

K  1 

calculations  are  implemented  by  distributed  processing,  the  numerators 
of  Equation  (2-5),  for  k  **  1,2,. ...N,  can  be  calculated  in  parallel, 
along  with  each  of  the  N  separate  filters.  An  estimate  of  <i  can  be 
calculated  but  is  not  required  to  calculate  state  estimates.  The 
conditional  covariance  of  x(t^)  is 


P(t+) 


1=1  +  Uc(tt} 

•  ■  ^(ti}3  T} 


-  x(<)  1 


(2-7) 


where  P^Ct^)  is  the  state  error  covariance  calculated  by  the  Kalman 
filter  based  on  a^.  If  P^t*)  is  desired,  it  must  be  computed  on-line 
because  it  depends  on  the  actual  measurement  history;  however,  it  is  not 
essential  to  the  operation  of  the  on-line  algorithm.  As  seen  in 
Equation  (2-6),  the  state  estimate  resulting  from  this  development  is 
the  probabilistically  weighted  average  of  the  state  estimates  generated 
by  N  separate  Kalman  filters. 


2.2  The  Multiple  Model  Adaptive  Filter 


The  multiple  model  adaptive  filtering  (MMAF)  algorithm  is  based  on 
the  previous  development.  As  depicted  in  Figure  2-1,  the  MMAF  is 
composed  of  N  separate  Kalman  filters,  each  based  on  a  discrete  value, 
a^,  of  the  parameter  vector.  At  time  t^,  the  measurement  is 
processed,  generating  residuals  ,  ^(t^)  *  *  *  *rN^fci^  ’  res^ua^s 
are  used  to  evaluate  the  hypothesis  conditional  probabilities,  p^t^) , 
according  to  Equation  (2-5).  The  probability  p^(t^_^)  needed  for  this 
iteration  is  maintained  from  the  last  sample  period  and 


(2„)"/2  [Afc{ti) | 1/2 


exp  {*} 


{*)  -  (-1/2 


(2-8) 


where  m  is  the  number  of  measurements  and  A^(t^)  is  available  from  the 


kth  filter  as 


W  ■  *  4V 


(2-9) 


Finally,  the  adaptive  state  estimate  is  generated  as  the  probablis- 
tically  weighted  average  seen  in  Equation  (2-6) ,  portrayed  as  the  output 
of  the  summing  junction  in  Figure  2-1, 

One  expects  that  the  filter  that  most  nearly  represents  the  true 
target  motion  will  produce  a  true  residual  covariance  very  close  to  its 
internally  computed  A^(t^) .  The  "mismatched"  filters  should  produce 
larger  residuals  than  anticipated,  so  that  Equation  (2-5)  leads  to  a 
heavier  weighting  of  the  estimates  produced  by  the  "best"  filter.  The 
algorithm's  performance  depends  on  the  existence  of  significant 
differences  between  the  residuals  of  the  "correct"  filter  and  the 
"mismatched"  filters.  These  significant  differences  can  be  enhanced  by 
specifically  tuning  each  filter  for  its  best  performance  against  a 
target  trajectory  that  matches  its  internal  dynamics  model.  The  common 
practice  of  "conservative"  tuning  of  a  single  nonadaptive  Kalman  filter 
by  adding  additional  dynamics  pseudo-noise  should  be  avoided,  since  it 
tends  to  blur  the  distinctions  between  the  estimates  and  residuals  based 
on  different  models. 

Additionally,  the  calculated  probabilities  should  have  an 
artificially  enforced  lower  bound.  By  bounding  the  lower  value,  the  p^ 
of  a  "bad"  filter  is  prevented  from  converging  to  zero,  which  would 
effectively  remove  that  filter  from  the  bank  and  prevent  the  MMAF  from 


responding  as  accurately  or  rapidly  to  future  changes  in  the  parameter 
values.  After  bounding,  the  p^  values  are  rescaled  so  that  the 
probabilities  still  sum  to  one.  The  selection  of  the  lower  bound  value 
can  significantly  affect  the  accuracy  of  the  estimate  due  to 
inappropriate  weighting  of  the  "mismatched"  filters.  A  bound  of  0.001 
is  selected  for  this  effort. 

In  this  application,  the  MMAF  consists  of  a  bank  of  either  two  or 
three  elemental  filters.  Each  elemental  filter  has  the  same  form; 
however,  they  are  tuned  for  different  assumed  target  trajectories.  Each 
filter  acquires  data  from  a  wide  tracking  window  (24-by-24  pixels)  or  a 
narrow  tracking  window  (8-by-8  pixels) ,  allowing  the  filter  to  maintain 
lock  on  a  highly  dynamic  target  or  providing  greater  resolution  for  more 
benign  target  motion,  respectively.  Specifically,  the  two-filter  MMAF's 
first  filter  is  optimized  for  a  benign  trajectory  and  an  associated 
small  tracking  window,  and  the  second  filter,  for  a  20  g  pull-up  target 
maneuver  with  a  larger  tracking  window  being  used.  The  three-filter 
MMAF  contains  the  first  two  elemental  filters  plus  a  narrow-FOV  filter 
tuned  to  perform  well  against  10  g  maneuvers.  By  changing  the 
probabilistic  weights,  the  MMAF  adaptively  changes  tracking  window  size 
and  dynamics  model.  The  current  two-filter  MMAF  has  a  very  coarsely 
discretized  parameter  space;  this  results  in  poorer  performance  against 
trajectories  which  are  not  modeled  by  one  of  the  elemental  filters.  A 
finer  discretization,  yielding  three  or  more  elemental  filters,  plus  the 
added  ability  to  restart  a  diverging  element  filter  (using  estimates 
from  the  nondiverging  filter(s)  to  accomplish  the  "restart")  should 
produce  more  accurate  tracking  than  purposefully  detuning  the  elemental 
filters  to  prevent  divergence  of  the  MMAF  when  tracking  a  20  g  maneuver. 
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Maximum  a  posteriori  (MAP)  state  estimation  is  based  on  calculating 

the  mode  of  f  ,  .  Under  the  standard  linear  system,  Gaussian  noise 

assumptions  for  a  Kalman  filter  development,  the  conditional  mean, 

x(t*),  is  the  mode.  However,  due  to  the  uncertain  character  of  the 

parameter  vector  one  must  consider  the  joint  conditional  density 

function,  f  In  this  case,  it  is  reasonable  to  assume  that  the 

x,a|Z 

conditional  mean  (or  mode) ,  calculated  by  the  filter  based  on  the  most 
probable  a^  (the  one  with  the  highest  p^)  is  a  good  approximation  lo  the 
true  mode  of  f  .  .  This  is  not  just  an  approximation  if  a_  is  truly 
discrete-valued . 

The  MAP  multiple  model  filter  calculates  p^(t^)  as  shown  in  Section 
2.2.  However,  unlike  the  Bayesian  MMAF,  the  MAP  estimate  is  not  the 
optimally  weighted  average  of  the  elemental  filter  estimates.  Instead, 
it  is  assumed  that  the  single  filter  with  highest  p^(t^)  is  “ost 

likely  representation  of  the  true  target  behavior.  The  estimate, 
x^(t^),  produced  by  this  single  filter  is  used  as  the  MAP  MMAF  estimate 
until  another  filter  becomes  the  "best"  match  to  the  actual  target 
dynamics . 

The  MAP  approach  to  estimation  has  positive  and  negative  aspects  as 


compared  to  the  Bayesian  approach.  It  is  expected  that  the  MAP  estimator 
will  respond  more  rapidly  to  changes  in  target  maneuver  characteristics. 
Moreover,  when  the  true  target  motion  is  well  matched  to  one  of  the 
filters  in  the  bank,  estimation  accuracy  is  not  degraded  significantly 
by  including  information  from  the  "mismatched"  filters.  However,  since 
the  need  to  reduce  computational  loading  requires  a  coarse 


discretization  of  the  parameter  space,  the  true  target  behavior  often  is 
not  well  matched  to  any  one  filter  in  the  bank.  In  this  case,  the 
estimator  operates  with  an  inappropriate  dynamics  model  and  can  be 
expected  to  provide  degraded  performance.  The  Bayesian  estimator 
provides  higher  accuracy  when  none  of  the  models  match  the  true  target 
behavior  because  of  its  probabilistically  weighted  averaging.  Each  of 
the  approaches  are  based  on  sound  reasoning,  leaving  experimental 
results  to  determine  which  approach  provides  the  best  results  for  a 
given  problem. 


III.  TRUTH  MODEL 


3. 1  Introduction 

The  truth  model  provides  the  simulation’s  standard  of  performance; 
the  standard  to  which  the  filter  estimates  are  compared  and  the 
"reality"  from  which  the  measurement  input  is  constructed.  As  such,  it 
should  provide  a  close  approximation  to  the  real  world  environment, 
particularly  in  the  areas  critical  to  the  problem  of  concern.  The 
processes  of  importance  in  this  investigation  are  atmospheric  jitter, 
target  dynamics  and  shape  effects,  and  background  and  FLIR  noises. 
Vibrational  effects  can  also  be  significant,  but  are  neglected  here, 
based  on  the  assumption  of  a  ground  based  weapon  system.  These 
processes  are  important  because  they  contribute  to  apparent  target 
motion  as  observed  by  the  sensor. 

Apparent  target  motion  is  described  using  an  x-y  coordinate  system 
in  the  FLIR  image  plane,  expressed  in  units  of  pixels  (a  pixel  is  20 
p  rads-by-20  prads) .  The  coordinates  of  the  observed  target  centroid 
are  the  sum  of  the  offsets  from  the  center  of  the  FOV  due  to  actual 
target  motion  and  the  offsets  due  to  atmospheric  disturbances 

XC  "  *D  +  XA 

where 

x  -  observed  centroid  x  coordinate 

\j 

Xp  -  x  coordinate  of  the  offset  due  to  dynamics 

x^  -  x  coordinate  of  the  offset  due  to  atmospherics 


and  similarly  for  y. 


Target  motion  in  the  truth  model  is  described  by  a  linear 


stochastic  differential  equation  [5:163] 


x(t)  -  F  x(t)  +  B  u(t)  +  G  w(t) 


where  the  state  vector,  x(t) ,  is  composed  of  the  target  position  and 
atmospheric  states.  The  Bki  term  consists  of  deterministic  velocities 
applied  to  move  the  target  in  a  preselected  inertial  trajectory.  The 
white  Gaussian  noise  term  drives  the  atmospheric  jitter  "shaping 
filter"  to  generate  xA  and 

The  measurements  used  by  the  filter  algorithm  are  derived  from  the 
truth  states  using  the  target  image  intensity  function  developed  later 
in  this  chapter.  Values  derived  from  this  intensity  function  are 
corrupted  by  adding  spatially  correlated  and  temporarily  uncorrelated 
noise,  accounting  for  FLIR  and  background  noises.  The  64-element 
vector  thus  constructed  is  compared  to  a  template,  derived  from 
earlier  data,  by  the  correlator  portion  of  the  tracking  algorithm  to 
produce  the  offset  values  used  as  measurement  input  to  the  linear 
Kalman  filter.  The  measurement  is  related  to  the  state  variables  by 
the  following  discrete-time  linear  stochastic  equation 


zC^)  »  H  x(t±)  +  v(t1) 

The  target  model  and  the  measurement  model  incorporated  in  the 
same  simulation  truth  model  are  described  in  this  chapter,  which  covers 


target  dynamics,  atmospheric  effects,  the  various  coordinate  frames, 


the  transformation  of  inertial  motion  to  the  FLIR  image  plane 
coordinates,  and  the  inertial  target  trajectories  used  for  this  effort. 
The  target  intensity  function,  the  target  image  projection,  and  noise 
effects  are  defined  in  greater  detail  as  well.  Significant  aspects  of 
the  truth  model  could  be  replaced  with  actual  data;  however,  this 
approach  of  simulating  a  truth  model  provides  control  over  many  aspects 
not  typically  available  from  test  data. 

3.2  Truth  State  Model 

The  truth  model  describes  the  apparent  motion  due  to  dynamics  and 
atmospheric  effects.  The  target  dynamics  model,  developed  by  Harnly 
and  Jensen  [3],  provides  the  true  location  of  the  target  center  of  mass 
in  the  FLIR  image  plane  horizontal  and  vertical  directions. 
Additionally,  this  section  presents  the  equations  which  translate  a 
simulated  target  inertial  trajectory  into  two  dimensional  motion  on  the 
FLIR  image  plane. 

In  order  to  evaluate  specific  trajectories,  a  deterministic  model 
provides  the  target  location  time  history.  The  first  two  elements  of 
the  state  vector,  and  are  the  dynamic  states  in  the  horizontal 
and  vertical  directions  in  the  FLIR  image  plane.  Because  of  the  large 
distances,  the  azimuth  angle,  a  ,  and  elevation  angle,  3  ,  are 
essentially  the  linear  coordinates,  x  and  y,  of  the  target  centroid 
position  in  this  FLIR  image  plane.  More  details  regarding  this 
relation  are  provided  later  in  this  chapter.  In  order  to  conform  to 
the  linear  differential  equation  format,  the  deterministic  inputs  must 
be  a  and  3  instead  of  a  and  8  directly.  Based  on 


the  model  s  dynamic  motion  can  be  expressed  as 


4(0  =  u(t)  =  [ a(t)  6  (t)]T  (3-1) 

The  atmospheric  disturbances,  developed  by  Mercier  [6],  are 
modeled  as  third  order  Gauss-Markov  processes  in  both  the  horizontal 
(x)  and  vertical  (y)  image  plane  directions.  In  the  x-direction  this 
yields 


where 

■  unit  strength,  zero-mean,  white  Gaussian  noise 
K  =  system  gain,  adjust  for  desired  atmospheric  RMS  value 
A  =  break  frequency;  14.14  rads/sec 
B  =  break  frequency;  659.5  rads/sec 
xA  *  output  of  the  shaping  filter 

Thus,  the  atmospheric  jitter  is  represented  in  x  and  y  directions 
by  the  stochastic  differential  equation 


4(t)  =  4  — A(t)  +  4  ^(t) 


(3-2) 


x^(t)  =  six  atmospheric  noise  states 
=  atmospheric  plant  matrix 
=  atmospheric  noise  input  matrix 
w^Ct)  =  vector  of  white  Gaussian  noise  inputs  with  statistics 


E{w^(t) }  =  0 

E(wA(t)wAT(t+ t)  }  =  (^(t)  <$(x) 


The  atmospheric  states  are  augmented  with  the  dynamic  states,  producing 
an  augmented  state  differential  equation  with  components  given  by 
Equations  (3-1)  and  (3-2) .  The  augmented  system  is  then  converted  to 
the  equivalent  discrete  time  form  [3].  The  propagation  of  the  target 
motion  takes  the  form 
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where 

x(t.)  =  state  vector  (two  dynamic  states  and  six  atmospheric 

1  states) 

B, (t  )  =  2-by-2  dynamics  input  matrix  a  I  fit,  where 

1  “-'i.i-'i 

Ud(ti)  =  piecewise  constant  function  (between  sample  times) 

evaluated  at  midpoint  to  approximate  the  integral 

of  a(t)  and  S(t)  from  t.  to  t.  . 
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The  state  transition  matrix  is 


represents  the  Cholesky  square  root  of  Q 
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Equation  (3-3)  is  developed  in  detail  in  reference  [3], 

As  stated  earlier,  azimuth  velocity,  a  (t) ,  and  elevation  velocity, 
3(t),  are  derived  from  deterministic  inertial  velocities.  The  inertial 
velocities  must  be  projected  into  the  FLIR  image  plane  coordinates. 


x^,  z_^.  =  inertial  axes 

r  =  range  from  tracker  origin  to  target 

r,  =  horizontal  range 
n 

v  =  target  inertial  velocity 
a  =  azimuth  angular  displacement 
6  =  elevation  angular  displacement 

The  geometry  associated  with  azimuth  direction  is  shown  in  Figure  3-2. 


Figure  3-2.  Azimuth  Geometry 

from  Figure  3-2 

a  (t)  =  tan  *  z^(t)  (rads)  (3-4) 

xx(t) 

L  _ 

and  so 

a(t)  =  xI(t)z].  -  zI(t)xI(t)  (rads/sec)  (3-5) 

zx2(t)  +  xT^(t) 

The  azimuth  velocity  from  Equation  (3-5)  is  in  rads/sec,  which  must  be 

_6 

converted  to  pixels/sec  by  dividing  by  20  x  10  rads/pixel  [3:33]. 


Similarly,  Figure  3-3  illustrates  the  geometry  involved  in 


a 

\'s 


computing  the  elevation  velocity. 


Figure  3-3.  Elevation  Geometry 


where 


r(t)  =  range  »  [x^t)  +  y-j-2 ( t)  +  z^Ct)  l1^2 

2  2  1/2 

r.  (t)  =  horizontal  range  =  [Xj  (t)  +  Z.J.  (t)  ] 


and 


3 (t)  *  tan 


-1  r 


yx(t) 


rh(t) 


(rads) 


so  that 


r2(t) 


where,  from  the  r,  (t)  expression  above,  it  can  be  seen  that 
h 


r^(t)  =  xI(t)xI(t)  +  z1(t)Zj'(t) 


rh(t) 


Once  again,  the  velocity  must  be  converted  to  pixels/sec. 


(3-6) 


3  (t)  =  rli(t)yI(t)  -  y1(t)r‘h(t)  (rads/sec)  (3-7) 


By  substituting  Equations  (3-5)  and  (3-7)  into  Equation  (3-1),  the 


truth  model  generates  target  motion  for  the  desired  trajectory. 


3.3  Target  Trajectories 

A  number  of  deterministic  trajectories  are  available  which 
incorporate  several  different  maneuver  options.  These  trajectories  are 
designed  to  provide  realistic  target  behavior  with  fairly  simple 
models.  The  basic  equations  are  described  in  detail  by  Millner  [8]. 

Trajectory  one  -  This  trajectory,  as  shown  in  Figure  3-4,  is  very 
benign.  The  target  flies  a  constant-heading,  straight-and-level 
course.  The  inertial  velocity,  v^.,  is  constant  for  the  entire  maneuver 
and  is  parallel  to  the  x^-z^  plane. 


Figure  3-4.  Trajectory  One 

Either  wings-level  flight  or  a  constant  roll-rate  maneuver 
(positive  roll  rate  defined  by  the  right  hand  rule)  may  be  simulated. 

Trajectory  Two  -  In  order  to  evaluate  the  filter  response  to  more 
dynamic  behavior,  this  trajectory  simulates  a  constant-g  pull-up.  The 
target  initially  flies  trajectory  one  (wings  level,  one-g  flight)  until 
t  =  2.0  seconds,  allowing  the  filter  to  obtain  good  position  estimates 


before  the  maneuver  is  initiated. 


I 


Figure  3-5.  Trajectory  Two 

The  pull-up  is  started  with  a  step  change  to  the  pitch  rate. 
While  this  is  an  unrealistic  behavior,  it  is  more  harsh  than  maneuvers 
encountered  in  the  real  world,  and  the  tracker  should  perform  at  least 
as  well  against  real  world  targets. 

Trajectory  Three  -  This  trajectory  contains  two  maneuver  changes, 
providing  a  means  to  evaluate  tracking  performance  against  a  target 
that  begins  and  ends  a  pull-up  maneuver.  As  with  trajectory  two,  a 
constant-g  pull-up  is  executed;  however,  the  pull-up  is  terminated  at 
t  «  3.5  seconds,  prior  to  the  end  of  the  simulation.  The  target 
acceleration  is  impulsively  set  to  zero;  the  inertial  velocity  the 
target  has  at  that  point  remains  constant  until  the  end  of  the 
simulation. 

Trajectory  Four  -  This  trajectory  displays  motion  in  all  three 
inertial  directions.  Like  trajectory  two,  the  target  initially  flies 
straight  and  level  and  then,  at  t  =2.0  seconds,  the  target  performs  a 
constant-g  pull-up,  but  in  the  -z  ^  direction  rather  than  the  +y^ 
direction  (see  Figure  3-5).  The  projected  target  image  on  the  FL1R 
image  plane  changes  more  dramatically  in  this  out-of-plane  maneuver. 
This  variation  allows  investigation  of  the  filter's  response  to  dynamic 
changes  in  the  orientation  and  spacing  of  the  individual  hot-spots  that 


compose  the  target  intensity  profile  on  the  FLIR  image  plane. 
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3.4  Coordinate  Frames 

The  dynamics  of  various  parts  of  the  target  relative  to  its  center 

of  mass  must  be  known  to  project  the  target  image  accurately  onto  the 

FL1R  plane.  This  projection  is  required  to  calculate  the  intensity 

values  for  measurement  inputs.  Definition  of  two  additional  coordinate 

frames  facilitates  these  calculations. 

Target  Frame  -  The  origin  of  the  target  frame  is  the  target's 

center  of  mass.  The  first  axis  is  coincident  with  the  velocity  vector. 

Perpendicular  to  the  first,  the  second  axis  points  out  the  right  side 

of  the  target.  The  final  axis  completes  the  right  hand  system  and 

points  out  the  underside  of  the  target  fuselage.  In  this  simulation, 

the  multiple  hot  spots  are  assumed  to  lie  in  the  e^-e^  plane.  The 

target  coordinate  system  is  expressed  as  unit  vectors  e  ,  e  ,  and  e 

~v  pv  — ppv 

(v:  direction  of  the  velocity  vector,  pv:  perpendicular  to  velocity, 
ppv:  perpendicular  to  both). 

a  -  8  plane  -  The  origin  of  this  frame  is  also  the  target  center  of 
mass.  One  basis  vector,  e  ,  is  aligned  with  the  true  line  of  sight 

T 

from  the  tracker  (located  at  the  inertial  coordinate  system  origin). 
The  plane  of  interest  is  defined  by  unit  vectors  e_a  and  e^  which  are 
rotated  from  the  inertial  frame  by  the  angles  a  and  8  ,  as  shown  in 
Figure  3-1. 
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3.5  Measurement  Model 


The  measurements  developed  for  the  tracking  filter  are  derived 
from  the  projection  of  the  target  intensity  function  onto  the  FLIR 
image  plane.  FLIR  and  background  noises  corrupt  the  intensity 
function.  At  close  range,  targets  can  be  well  modeled  as  the  sum  of 
bivariate  Gaussian  functions  with  elliptical  contours  [3].  As  shown  in 
Figure  3-6,  the  apparent  target  image  eventually  resolves  into  separate 
hot  spots. 

CENTROID  OF 
APPARENT  TARGET 
EQUAL-INTENSITY  INTENSITY  PROFILE 


Figure  3-6.  Apparent  Target  Intensity  Pattern  on  FLIR  Image  Plane 


This  simulation  used  three  identically  distributed  hot  spots,  each  of 
which  is  described  by  the  following  intensity  function. 


1  [x'y’xpeak(t)’ypeak(t)]  *  lmax  eXp  {“°-5  [x-Xpeak}  ^peak* 1 

'  It''1  l-W  <y-ypeak>IT 

(3-8) 


where 


I  -  maximum  intensity  of  the  hot  spot 
max 

t  ,  y  -  coordinates  of  the  peak  intensity  of  the  hot  spot 

p63K  p6dK 

■>  2 

-  matrix  whose  eigenvalues  are  a  and  a  ,  which 
are  the  dispersions  ofV  the  pv  elliptical 
constant-intensity  contours  in  the  target  frame, 
and  whose  eigenvectors  define  the  orientation  of 
the  ellipse  principal  axes 


The  x-  and  y-coordinates  in  this  function  are  calculated  in  pixels 
relative  to  the  center  of  the  tracker  field  of  view. 

The  location  of  the  hot  spots  are  expressed  in  the  target 
coordinate  frame.  The  intensity  function  centroid  for  single  hot  spot 
targets  are  assumed  coincident  with  the  target  center  of  mass.  In  this 
study,  the  multiple  hot  spots  are  distributed  as  shown  in  Figure  3-7 


Figure  3-7.  Hot  Spot  Distribution 


Hot  spot  one  is  located  at  1.0  meters  in  the  +e  direction;  two  and 


three  are  located  at  +0.5  meters  in  the  e  direction.  It  is  assumed 

-  -pv 

that  the  semi-major  axes  of  the  ellipses  are  parallel  and  aligned  with 
the  velocity  vector.  The  velocity  vector  remains  out  the  nose  of  the 
target  for  one  entire  simulation.  While  these  are  not  realistic 
assumptions,  they  facilitate  the  simulation  of  target  dynamics,  and  the 
tracking  algorithm  can  be  expected  to  perform  as  well  against  live 
targets . 

The  measurements  for  both  single  and  multiple  hot  spot  cases  are 
the  average  intensity  values  from  each  pixel  of  an  8-by-8  tracking 
window.  This  intensity  is  the  sum  of  each  hot  spot  contribution  and 
background  and  FLIR  noises.  The  measurement  value  for  pixel  kl: 


^  fl/AP  j1  m[x*y*xpeakm(ti)’ypeakm(ti)ldxdy}  +  vkl(ti> 

pixel  kl  (3-9) 


where 

I^t*]  =  intensity  function  of  the  m^  hot  spot  of  M  total  hot  spots 

z^(t^)  =  output  of  the  klth  pixel  (kth  row,  1th  column)  at  time  t^; 

the  average  intensity  at  that  pixel  as  sensed  by  a  detector 
in  the  FLIR  image  plane. 


area  of  one  pixel 


,  th 


(x. 


(x,y)  ~  coordinates  of  any  point  within  the  kl  pixel 

.  ,y  ,  )  =  location  of  one  peak  of  the  m*”*1  intensity  function 
peakm  'peakm  r 

at  L  , 


v^(t^)  =  additive  FLIR  and  background  noises  for  the  kl*'*1  pixel 


3.6  Target  Image 

During  the  course  of  the  simulation,  the  target  size  and  hot  spot 
distribution  on  the  target  will  not  change;  however,  as  the  target 
maneuvers  and  approaches  the  sensor,  the  target  image  on  the  FLIR 
sensor  does  vary.  In  order  to  simulate  these  variations,  the  target 
image  at  any  given  time  is  referenced  to  a  previously  defined  standard 
image.  The  target  image  is  developed  in  the  target  frame  and  then 
projected  onto  the  a  -  8  plane,  described  in  Section  3.3;  the  geometry 
is  illustrated  in  Figure  3-8. 


Figure  3-8.  Image  Projection  [9:11-23] 

The  reference  image  is  defined  with  the  target  flat  in  a  plane 
perpendicular  to  the  tracker  line  of  sight.  Any  other  orientation 
would  produce  a  smaller  image.  The  following  expressions  relate  the 
current  image  size  to  the  reference  image  based  on  the  current  range 


and  velocity  of  the  target  [9:11-24]. 


where 


i 


a  =  a  (p  /p) 

pv  pvo  o 

°v  51  (po/p)  t?pvo  +  <pvo  ■  V  COS  £  1 
-  0nv  {1+  P— 1L0S)  /— I-I 


a  vo’  ^>vo  =  the  dispersion  of  the  reference  target  image 

along  the  major  and  minor  axes  of  the  radiation 
ellipsoid,  parallel  and  perpendicular  to  the 
velocity  vector. 

av’°pv  *  The  current  dispersions  of  the  target  image 

p  =  reference  range  from  the  sensor  to  the  target 
o 

p  =  current  range 
v  =  inertial  velocity 

(v^L0S^  =  projection  of  onto  the  a  - B plane 

£  =  angle  between  the  inertial  velocity  vector  and 
the  o-B  plane,  shown  in  Figure  3-8 

A.R  a  (j  Jo 

vo  pvo  -  maximum  aspect  ratio  of  the  reference 
image 


The  image  resulting  from  the  above  equation  is  defined  in  the 
target  frame  coordinates.  Both  the  coordinates  of  the  hot  spot  centers 
and  the  dispersion  matrix  must  be  transformed  into  a  -  6  plane 
coordinates  [9:11-25].  From  Figure  3-8 

cos  9  *  o  (t)/[v^LQS3 
sine  -e  <t)/[viL0S] 


where  (v.  „)  is  the  magnitude  of  the  velocity  perpendicular  to  the 

XLUu 

.  2  •  2  1/2 

tracker  line  of  sight,  defined  as  [ a(t)  +  8(t)  ]  .  The  hot  spot 

coordinates  are  transformed  to  a  -  6  coordinates  by 


••  c-  .s  .•  .  - , 


X 


cos 9  -sin0 


X 


in0  cos0 


-  J  target  frame  =  A  x 


The  dispersion  matrix  is  transformed  using 


—  a6  =  ^  P  A* 


The  inverse  of  IP  is  required  for  Equation  (3-8) .  This  is 
conveniently  obtained  by  using  the  equivalent  expression 


P  1  =  A(P_1)AT 

a8 

These  values  are  used  in  Equation  (3-8)  to  calculate  the  intensity 
values  for  the  measurement  update. 


3.7  Spatially  Correlated  Background  Noise 

The  noise  term,  v^(t^)  ,  in  Equation  (3-9)  contains  spatially 

correlated  background  noise  [3].  The  correlation  distance  is  about  two 
pixels;  the  noise  is  modelled  by  non-zero  circularly  symmetric 
correlation  between  each  pixel  and  its  two  closest  neighbors  in  all 
directions. 

The  measurements  are  arrayed  as  a  64  element  vector  (64  pixels  in 
the  tracking  window) .  The  additive  noise  term  is  modeled  as 


v(t±)  =  ^v'(ti) 


where  v.'  (t^)  is  a  vector  of  64  independent,  discrete-time  zero-mean, 


white  Gaussian  noise  processes  with  a  variance  of  one.  The  strength  of 
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the  resulting  noise  process,  v^(t^) ,  is  E{v(t^)\^  (t..)}  =  R5„.  This 
64-by-64  matrix  describes  the  spatial  correlation  between  pixels  and  is 
presented  in  detail  by  Harnly  and  Jensen  [3]  and  Kozemchak  [4],  The 
correlation  terms  not  involving  first  or  second  neighbors  of  a  given 
pixel  are  essentially  zero.  The  noise  term  added  to  the  average 
intensity  value  of  each  pixel,  forms  the  simulation  measurement  data. 
The  adaptation  for  a  wide-FOV  is  discussed  in  Section  4.4. 

3.8  Summary 

This  chapter  presents  the  truth  model  state  propagation.  Also 
included  are  the  transformation  from  the  deterministic  inertial  target 
trajectories  to  target  motion  on  the  FLIR  image  plane.  Finally,  the 
measurement  model  is  developed.  The  next  chapter  describes  the  linear 
Kalman  filter/correlator  tracking  algorithm  and  its  use  in  a  multiple 


model  filter. 


IV.  TRACKING  ALGORITHM 


4. 1  Introduction 

The  linear  Kalman  filter/correlator  tracking  algorithm  that  is 
the  elemental  unit  in  the  multiple  model  adaptive  filter  (MMAF)  was 
developed  by  Rogers  [II]  and  Millner  [7],  It  uses  a  linear  Kalmar 
filter  to  provide  estimates  of  target  position,  velocity,  and 
acceleration  and  atmospheric  disturbances.  An  enhanced  correlation 
algorithm,  using  the  estimated  target  shape  as  a  template,  provides 
the  Kalman  filter  with  pseudo-measurements.  These  pseudo¬ 
measurements  are  offset  distances  from  the  center  of  the  FLIR  FOV, 
modeled  as  a  linear  combination  of  the  filter  state  variables.  The 
linear  formulation  allows  many  terms  of  the  Kalman  filter  to  be 
precomputed,  thus  reducing  on-line  computational  loading. 


The  equations 

necessary 

to  propagate  and 

update  the 

state 

estimates 

are 

presented  in 

this  chapter. 

Additionally, 

the 

processing 

of 

the 

template 

and  the  correlation  algorithm 

are 

described.  The  linear  Kalman  filter/correlator  is  the  basic  element 
in  a  bank  of  similar  filters  used  by  the  MMAF,  as  discussed  in 
Chapter  II.  The  chapter  concludes  with  the  specific  variations  in 
each  elemental  filter  used  to  define  the  MMAF. 

4 . 2  State  Space  Model 

The  linear  Kalman  filter  models  target  acceleration  and 
atmospheric  jitter  position  as  stationary,  first-order  Gauss-Markov 


processes.  In  the  truth  model,  described  in  Chapter  III,  the 
atmospheric  jitter  position  was  represented  by  a  third  order  model; 
however,  the  high-frequency  double  pole  has  little  effect  on  the  low 
frequency  characteristics  of  the  jitter  and  is  neglected  in  the 
filter  model  [5] . 


The  filter  states  are  the  target  position  (x^,yD),  velocity 

(v  ,v  )  ,  and  acceleration  (a  ,a  ),  and  iitter  position  (x.  ,y.),  in 
xy  x  y  A  A 

each  FLIR  image  plane  direction. 


5*  =  [XD  yD  Vx  ax  ay  XA 


(4-1) 


The  state  equations  are 


where 


=  correlation  time  for  target  acceleration 
x^p  =  correlation  time  for  atmospheric  jitter 


zero-mean  white  Gaussian  noise  processes  whose 
strengths  depend  on  tuning  results 
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Identical,  independent  models  are  used  to  represent  the  effects  in 
the  x-  and  y-directions  of  the  FLIR  image  plane. 

The  above  relationships  can  be  written  as  a  state  vector 
differential  equation  in  standard  form 


where 


*F(t)  =  FpX/t)  +  Gf  w^Ct) 


(4-2) 


Fy  =  The  time  invariant  system  plant  matrix  which  is 


'  DF 

0Ut  -1/t , 


the  time  invariant  system  noise  output  matrix 


wF(t)  =  twDx'  WDy’  WAx*  WAy^  ’  noise  vector  of  mutually 

independent  zero-mean  white  Gaussian  noise 
processes 


KfWjXtJVy  (t+r)}  *  Qp  >5(t) 
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ctdf  =  target  acceleration  variance  and  ms  value 
2 

o.-p  =  atmospheric  jitter  position  offset  variance  and  ms 
value 

The  equivalent  discrete  time  equation  [5:173]  is 


where 

XyCt^)  *  filter  state  vector  at  time,  t^ 

WFp(tj)  *  discrete-time  zero-mean  white  Gaussian  noise  of 
covariance  *  I 


*F(ti+P  =  *F(ti+l  "  ti)^F(ti)  + 


The  state  transition  matrix  is  [13]: 


$  (t 
Fv  l+l 


1  0 
0  1 
0  0 
0  0 
0  0 
0  0 
0  0 
0  0 


t  0 
0  t 
1  0 
0  1 
0  0 
0  0 
0  0 
0  0 


J1  0 

0  J1 

J2  0 

0  J2 

J3  0 

0  J3 

0  0 

0  0 


0  0 

0  0 

0  0 

0  0 

0  0 

0  0 

J4  0 

0  J4 


where 

J1  ■tDF14t-TDF(1-e!tP<-4t/'DF»l 
J2  -  tDF(l-exp(-At/TDr)] 

J3  *  exp(-At/xDF) 

J4  =  exp(-At/xAf,) 
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The  growth  in  uncertainty  due  to  dynamic  driving  noise  is  [8] : 
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where 
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(2  TDFAt3/3]  -  [2  tdf2  At2]  -  [4  tdf3 
[2  tdf3  At]  -  [Tdf4  exp (-2  At/x^)  ]  ♦  x 


At 

4 

DF 


exp(-At/TDF) ] 


^13  “  °DF  ^TDF  At  1  +  f2  TDp  At  exp(-At/xDF) ]  +  [xDp  ] 

-  t2  *DF3  exp(-At/xDF)]  -  [2  xDF2  At]  +  [x^3  exp(-2  At/x^)  ] 

Ql5  =  aDF  t“2  x^F  At  exp(-At/XpF) ]  +  [x^F  ]  -  [x^F  exp(-2At/TpF) ] 

Q33  =  aDF2  12  tDF  Atl  *  [3  tDF2]  +  [4  tDF2  exp(-At/xDF) ] 

-  [xDF2  exp (-2  At/xDp)] 

Q35  =  CTDF2{TDF  '  12  tDF  exp(-At/TDF)]  +  [xDp  exp(-2At/xDF) ]  > 

Q55  =  aDF2  { 1  ’  exp("2  At/TDF)} 

Q77  =  (  1  ~  exp (-2  At/x^)} 


Using  the  terms  defined  above,  the  state  estimate  and  covariance 
are  propagated  as  follows: 


<tF(ti+l,ti)  -^Ci^ 

(4-4) 

KV?  ‘ 

*F(ti+l»ti)i(tt) $F  (ti+l,ti)  +  SfD 

(4-5) 

(4-5) 


x(t^) 

£(ti+p 

p(t^) 


=  the  state  estimate  after  update  at  time  t^ 

=  the  state  estimate  prior  to  update  at  time  t^+^ 

=  the  conditional  covariance  matrix  after  update  at  time  t^ 

=  the  conditional  covariance  matrix  prior  to  update  at 


time  t 


i+1 


4.3  Measurement  Model 

The  Kalman  filter  update  process  uses  pseudo-measurements 
produced  by  an  enhanced  correlation  algorithm.  The  template  used  by 
the  correlator  is  developed  by  preprocessing  the  target  image.  The 
following  subsections  describe  the  target  image  processing,  the 
correlation  algorithm,  and  the  Kalman  filter  update  equations. 

4.3.1  Two-Dimensional  Fourier  Transform.  Many  of  the 

operations  required  to  perform  target  image  estimation  are  more 

easily  accomplished  in  the  Fourier  domain.  Additionally,  Fourier 

transforms  are  readily  accomplished  optically;  optical  implementation 

of  the  tracking  algorithm  would  greatly  reduce  the  required  computer 

resources.  The  Fourier  transform  of  a  complex-valued  function  of  two 

independent  variables,  g(x,y) ,  is  a  decomposition  of  g(x,y)  into  a 

linear  combination  of  functions  of  the  form  exp[j2ir(f  x  +  f  y)  ] .  The 

x  y 

transform  is  defined  as  a  double  integral  with  respect  to  the  spatial 
variables  x  and  y,  where  f^  and  f^  are  spatial  frequencies  [12]. 

In  this  case,  the  discrete  Fourier  transform  (DFT)  is  used 
because  the  FLIR  measurements  are  discrete  values  of  average 
intensity.  The  DFT  and  its  inverse  are: 
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N-l  N-l 

l  E  h(x,y)  exp[-j2  tt  (f  x  +  f  y)  ] 
x=0  y=0  x  y 


h(x,y) 


1/N 


N-l  N-l 

E  E  H(f  ,f  )  exp[+i2  it  (f  x  +  f  y)  ] 
f  f  x  y  r  J  x  yJ 

I _ *• _ n 


where 


H(f  f  ) 
x  y' 


transformed  function  in  the  spatial  frequency 
domain 


h(x,y)  =  function  in  the  spatial  domain 


f 


f 

y 


spatial  frequencies 


x,y  =  spatial  variables 

N  =  the  period  of  the  assumed  recurring  sequence  in 
both  directions;  thus  the  sequence  of  intensity 
values  is  discretized  into  an  N  x  N  pixel  array. 


The  8-by-8  tracking  window  is  padded  with  additional  data  to 
form  a  24-by-24  array  for  processing.  The  original  array  is  padded 
to  reduce  edge  effects,  aliasing,  and  leakage  conditions.  Data, 
available  from  the  larger  field  of  FLIR  data,  is  used  for  padding 
instead  of  zeros  to  prevent  artificial  edge  effects  [4:19].  The 
"large"  FOV  data  is  padded  with  zero  for  reasons  discussed  in  Section 


4.4. 


4.3.2  Fourier  Transform  Shifting  Property.  In  order  to  average 
the  target  image  temporally,  the  intensity  function  in  successive 
frames  of  data  must  be  centered  in  the  current  FOV.  The  intensity 
function  is  centered  on  the  filter's  estimated  target  centroid 
location  using  the  shifting  property  of  the  Fourier  transform  [9]. 
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The  Fourier  transform  shift  theorem  states  that  a  linear  phase 

shift  in  the  frequency  domain  corresponds  to  a  translation  in  the 

spatial  domain.  This  phase  shift  can  be  thought  of  as  a  cylindrical 

shift  due  to  the  assumed  periodic  nature  of  the  sampled  data.  Thus, 

using  this  property,  the  only  difference  between  the  centered  image 

and  a  translated  image  is  a  linear  phase  shift  proportional  to  the 

spatial  displacement  in  the  x-  and  y-directions.  If  F[g(x,y)]  = 

G(f  ,f  ),  then 
a  y 


F[g(x-a,y-b) ]  -  G(f  ,f  )  exp[-j2  (fa  *  f  b) ] 
x  y  x  y 


where 


a  =  spatial  translation  in  the  x-direction 
b  ■  spatial  translation  in  the  y-direction 


As  mentioned  earlier,  the  filter's  estimate  of  the  target  centroid 
location  relative  to  the  center  of  the  FOV  is  used  to  determine  the 
required  phase  shift  to  center  the  image  for  interframe  shooting. 

4.3.3  Exponential  Smoothing.  The  available  target  image 
measurements  are  corrupted  by  FLIR  and  background  noises.  It  is 
assumed  that,  for  the  applicable  sample  rates,  these  noises  tend  to 
vary  more  rapidly  than  the  target  intensity  pattern  from  sample 
period  to  sample  period  [11]. 

This  assumption  is  exploited  by  use  of  an  exponential  smoothing 
algorithm  to  time-average  the  data,  thereby  increasing  the  signal  to 
noise  ratio.  Exponential  smoothing  approximates  a  true  finite  memory 
average  without  requiring  the  storage  of  a  number  of  previous  frames 
of  data  [6:33],  The  applicable  equation  is 
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where 


2.U)  ■  ot£(t)  +  ( 1-a  )y_(t-l) 


£(t)  *  current  average  value 

2.(t)  =  current  data  frame 
2.(t-l)  =  previous  average  data  frame 
a  =  smoothing  constant,  0<ai  1 


The  appropriate  value  of  smoothing  constant  varies  according  to  the 
dynamics  of  the  target  image  variations.  A  rapidly  changing  image 
requires  heavier  weighting  of  more  recent  data  and  thus  a  higher 
value  of  a  ;  likewise,  slower  variations  require  smaller  values. 
Based  on  previous  studies,  a  is  set  at  0.1  for  this  effort  [13]. 

The  following  steps  summarize  the  template  derivation  process. 


1.  The  Fourier  transform  of  the  raw  FLIR  data  is  calculated 

2.  The  appropriate  negating  phase  shift  is  applied  to  center  the 
image  based  on  the  estimated  target  centroid  location 

3.  Temporal  smoothing  of  the  centered  data  is  performed 

4.  The  intensity  shape  function  is  evaluated  at  the  estimated  state 
after  the  control  application.  For  a  single  elemental  filter 


— (ti+l) 


*A(Vl> 


because,  as  stated  in  Section  1.3,  control  is  applied  so  as  to 
zero  the  predicted  dynamic  states.  Superscript  c  denotes  after 
application  of  the  control.  This  is  not  the  correct  relation 
for  the  MMAF,  as  explained  in  Section  4.4. 


4.3.4  Fast  Fourier  Transform  (FFT)  Correlation  Algorithm.  The 
correlator  used  in  this  effort  is  considered  enhanced  because  it 
compares  the  raw  data  with  a  template,  rather  than  the  previous  frame 
of  data.  The  template  is  the  estimated  target  intensity  function 
located  at  the  best  estimate  of  the  centroid  offset,  ,x(t^  ).  The 
FFT  is  used  to  perform  the  cross  correlation  by 

F[&(x, y) ]  =  £(fxfy) 

FU(x,y)]  -  L(f  f  ) 
x  y 

[£(x,y)  *  Ux.y)]  »  C(f  f  )  •  L*(f  ,f  )  (4-6) 

x  y  x  y 

where 

[*<*  ,y)  *  Ux,y)]  »  cross  correlation  of  the  two  dimensional 
spatial  sequences  £(x,y)  and  _l(x,y) 

L^*(f  ,f  )  -  complex  conjugate  of  the  Fourier  transform 
X  y  of  the  sequence  Kx,y) 

The  cross  correlation,  R(x,y),  is  obtained  by  taking  the  inverse 
FFT,  or  IFFT,  of  Equation  (4-6). 

The  resulting  cross  correlation  is  subjected  to  a  thresholding 
process  to  remove  false  peaks.  Any  element  of  R_(x,y)  that  is  less 
than  the  preselected  fraction  of  the  maximum  correlation  is  set  to 
zero.  Following  this,  a  centroid  summation  is  used  to  locate  the 
center  of  mass.  It  is  assumed  that  the  centroid  of  the  thresholded 
correlation  function  is  a  good  approximation  of  the  peak  location 
[11:53];  it  also  avoids  the  ambiguities  that  local  peaks  often  cause 
peak-finding  algorithms.  The  calculated  centroid  calculation  is  the 


correlator's  estimate  of  the  offset  of  the  target  from  the  center  of 
the  data  frame.  These  offsets  are  the  "measurements"  supplied  to  the 
linear  Kalman  filter. 

4.3.5  The  Kalman  Filter  Update  Equations.  The  linear  Kalman 
filter  updates  the  state  estimate  based  on  the  offset  values  from  the 
correlator.  The  appropriate  measurement  equation  is 

“  Hy  XF(ti)  +  (ti>  (4-7) 


where 


z^(t^)  =  the  offset  estimate  in  x-  and  y-coordinates 

produced  by  the  correlator,  based  on  the  filter 
predicted  centroid  location  due  to  dynamics  and 
atmospherics. 


XAD(ti^l  +(XAC(ti  i 


yAD(ti^J  +lyAC(ti). 


vFl(ti> 

vF2(ti^ 


H  =  10000010 

lo  1  0  0  0  0  0  1 


VpCt^^)  *  noise  produced  by  the  correlation  algorithm 

with  statistics  that  were  shown  empirically  to 
be  [11] 


ErVyU^}  =  0 


where 


R_,  =  .00436  0  " 

[  0  .00598 


The  standard  Kalman  filter  update  equations  apply: 


K(ti>  =  P F(t“)  Pp(t*>  HpT  +  Rp]_1  (4-8) 

Xp(t^)  =  2ip(ti)  +K(ti)[^(ti)  -HFxF(ti)]  (4-9) 

l F(tJ)  -  V*?  -  K(t.)  Hj,  Pp(tJ  (4-10) 

where  all  the  terms  have  been  previously  defined. 

4.4  Multiple  Model  Implementation 

The  theoretical  basis  for  the  multiple  model  adaptive  filter  is 
developed  in  Chapter  II;  this  section  concentrates  on  the  variations 
in  the  elemental  linear  Kalman  filter /correlators  included  in  the 
MMAF  bank  of  estimators.  The  single  filter  tracking  algorithm 
discussed  so  far  can  potentially  handle  a  target  maneuvering  at  a 

level  between  +20  g's,  with  tuning  appropriate  to  that  maneuver.  By 

2 

changing  the  values  of  t  and  a  used  to  model  acceleration  in  the 

Ur  Ufc 

filter  model,  each  filter  in  the  bank  can  be  tuned  for  different 
tracking  scenarios.  In  order  to  establish  a  filter  (or  more  than 
one)  that  can  maintain  lock  during  exceptionally  harsh  maneuvers,  a 
larger  tracking  window  is  incorporated  along  with  appropriate  tuning. 
The  larger  24-by-24  FOV  is  still  processed  as  an  8-by-8  data  array  in 
an  effort  to  limit  computational  loading.  However,  each  "large" 
pixel  is  60  y rads-by-60  yrads  instead  of  20  yrads  on  a  side.  This  is 
accomplished  by  making  each  scalar  intensity  data  value  the  average 
of  3-by-3  pixel  array.  The  same  FOV  center  is  maintained  for  both 
the  narrow  and  wide  tracking  windows . 


The  larger  FOV  requires  additional  modifications  to  the  basic 
algorithm.  The  measurement  noise  for  the  large  FOV  is  assumed  to  be 
spatially  uncorrelated  because  the  distance  between  "large"  pixel 
centers  is  greater  than  the  two  "small"  pixels  that  compose  the 
spatial  correlation  distance.  Also,  the  larger  FOV  makes  it 
reasonable  to  assume  the  target  intensity  value  at  the  edge  of  the 
tracking  window  is  essentially  zero;  this  makes  it  appropriate  to  pad 
the  data  array  with  zeros  instead  of  noisy  data,  when  estimating  the 
target  intensity  function  [13: XV— 9 1 .  Since  the  correlator  output  is 
in  terms  of  pixels,  "large"  or  "small",  the  pseudo-measurements  for 
the  large  FOV  case  are  multiplied  by  three,  allowing  the  filter 
portion  of  the  algorithm  to  remain  unaffected  by  the  different  size 
tracking  windows. 

This  effort  evaluates  three  MMAFs;  the  first  has  two  filters  in 
the  bank  and  the  other  two  MMAFs  add  a  third  filter;  one  of  these  is 
a  Bayesian  MMAF  and  the  other  is  a  MAP  MMAF  formulation.  The  first 
elemental  filter  for  each  MMAF  is  tuned  for  a  benign  trajectory, 
using  a  narrow  FOV  for  tracking  accuracy.  The  other  end  of  the 
expected  maneuver  range  is  covered  by  a  wide-FOV  filter  tuned  to 
track  a  20  g  target.  In  order  to  improve  accuracy  when  the  target 
motion  is  not  well  modelled  by  the  first  two  filters,  the  small  FOV 
third  filter  is  tuned  for  10  g  maneuvers.  The  tuning  results  are 
presented  in  Table  IV- 1,  below. 
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Table  IV- 1.  Filter  Tuning  Values 


2  2  4 

xDF(sec)  aDF  (pixels  /sec  ) 


MMAF 

1 

Filter  1 
Filter  2 

3.5 

1.5 

1000 

2000 

MMAF 

2 

Filter  1 

1.5 

1000 

Filter  2 

1.5 

1500 

Filter  3 

1.5 

2000 

2 

The  high  value  of  0pF  for  filter  1  in  the  first  MMAF  is  due  to  the 
coarse  discretization  of  the  parameter  space  and  is  required  to 
prevent  MMAF  1  from  diverging  when  tracking  targets  perform  maneuvers 
at  20  g's. 

Severe  maneuvers  can  cause  significant  differences  between  the 
position  estimates  of  the  MMAF  and  the  small  FOV  filters.  These 
differences  may  become  large  enough  to  cause  a  complete  cycle  of  the 
cylindrical  shift,  leading  to  divergence  of  the  small  FOV  filter. 
Due  to  the  artificial  lower  bounding,  the  Bayesian  MMAF  estimate 
continues  to  include  these  value,  ultimately  leading  to  divergence  of 
the  MMAF.  Additionally,  if  a  small  FOV  filter  is  allowed  to  diverge 
during  harsh  target  maneuvering,  it  will  not  be  capable  of  responding 
swiftly  if  the  target  motion  again  matches  its  internal  dynamics 
model  at  a  later  time.  To  avoid  this  problem,  a  reacquisition 
process  is  initiated  whenever  the  shift  of  the  narrow  tracking  window 
exceeds  a  magnitude  of  3.0  pixels.  The  divergent  filter  states  are 
reset  to  a  combination  of  the  nondivergent  filters'  states  or,  for 


the  MAP  estimator,  the  MMAF  state  estimate. 


Likewise,  for  the 


Bayesian  MMAF,  the  covariance  matrix  is  reset  to  the  combined 
covariance  values  of  the  nondivergent  elements;  the  MAP  MMAF  uses  the 
current  multiple  model  values  for  covariance.  The  conditional 
probabilities  are  left  at  the  current  values. 

Earlier,  it  was  pointed  out  that  the  estimated  target  intensity 
function  is  evaluated  at  the  predicted  target  centroid  position, 
after  control  application,  to  form  the  correlation  template.  With 
the  MMAF,  the  control  no  longer  zeros  out  the  predicted  offset  of 
each  elemental  filter  due  to  dynamics;  instead  the  FL1R  is  reoriented 
to  the  predicted  location  of  the  target  centroid  due  to  target 
dynamics,  based  on  the  MMAF  estimate.  The  intensity  function  is  now 
evaluated  at 

^DFk(ti+l)  +^AFk(ti+l)]  "  ^DMMF^i+P 

for  the  k  filter,  where  and  *AFk^Ci+l^  are  the  k 

filter's  target  dynamic  state  and  atmospheric  jitter  state  estimates 
(predictions)  ,  respectively,  and  x^^j,(t^+^)  is  the  MMAF  dynamics 
state  estimate  (prediction)  based  on  the  Bayesian  or  MAP  criteria 
developed  in  Chapter  II. 

4.5  Summary 

This  chapter  described  the  details  of  a  tracking  algorithm  based 
on  the  multiple  model  adaptive  filter  developed  in  Chapter  II.  The 
first  section  presented  the  propagation  equations  for  the  basic 


linear  Kalman  filter.  Next,  the  intensity  function  estimation 


process  used  Co  develop  the  template  for  the  enhanced  correlator  was 
outlined.  The  update  equations  completed  the  description  of  the 
linear  Kalman  filter/correlator  which  forms  the  basic  element  for  the 
multiple  model  filter  bank.  The  final  section  detailed  the 
variations  of  this  basic  algorithm  as  implemented  in  the  three  MMAFs. 


V.  ALGORITHM  TEST  SET-UP 


i 
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5. 1  Introduction 

This  chapter  presents  the  evaluation  methods  for  the  three 
tracking  algorithms  investigated  during  the  course  of  this  research. 
The  first  section  covers  the  derivation  of  tracker  statistics.  The 
format  of  the  performance  plots,  which  are  the  major  evaluation  tool, 
is  discussed  next;  followed  by  an  outline  of  the  parameter  values 
used  in  the  various  test  cases.  The  final  section  summarizes  the 
three  filters  and  defines  the  test  cases  used  to  generate  the  tracker 
performance  statistics. 

5.2  Tracker  Statistics 

Monte  Carlo  analysis  techniques  were  employed  to  generate 
tracker  performance  statistics.  Earlier  studies  established  that  ten 
runs  allow  sufficient  convergence  of  sample  error  statistics  to  the 
true  error  statistics  that  would  be  observed  with  an  infinite  number 
of  runs  [7,13]. 

Tracking  performance  is  determined  by  a  filter's  error  in 
estimating  the  target's  true  position  due  to  dynamics.  Consistently 
low  mean  errors,  along  with  low  error  variances,  demonstrate  the 
filter's  ability  to  track  the  target  through  a  wide  range  of 
maneuvers . 

The  sample  mean  errors  of  the  filter  position  estimates  are 


calculated  using 


(5-1) 


where 


-  <l/»)  J_1  [x*^)  -  xdfa(tt)l 

■  <1/N)  “  ,  W'l> 

n=l 


Exjj(t  )  *  sample  re  an,  error  (averaged  over  N  simulations) 
in  x-dynamics  position  at  time  t^ 

x,  (t.)  ■  truth  model  x-dynamics  value  at  time  t.  for  run  n 
an  I  l 

x,,.  (t, )  *  multiple  model  filter  estimated  x-dynamics  value 
din  i  ,  c 

at  time  t.  for  run  n 
i 

e  .  (t.)  =  number  of  Monte  Carlo  runs 
xdn  i 


The  sample  variance  of  the  error  is  calculated  by 


N 

“xd  (tl>  ■  l'/'8'1)!  lml  Sd/(tl>  -  Ex/(t,)  (5-2) 


where  the  terms  are  defined  above.  The  same  equations  also  apply  to 
calculation  involving  the  y-dynamics  position.  The  errors  are 
measured  in  FLIR  image  plane  coordinates  and  represent  offsets  from 
the  center  of  the  sensor  FOV.  Error  is  expressed  in  units  of  pixels; 
each  pixel  is  20  urads  on  a  side. 

Mean  error  and  standard  deviation  are  also  time  averaged  in 
order  to  provide  a  scalar  figure  of  merit  for  tracking  accuracy. 
This  temporal  averaging  occurs  for  a  short  period  after  transients 
have  died  out  until  just  before  the  first  maneuver,  and  from  just 
after  maneuver  transients  have  declined  until  a  time  close  to  the  end 
of  the  run.  Time  averaging  allows  presentation  of  data  in  a  compact 
tabular  form;  however,  these  figures  should  be  viewed  in  context. 


Misleading  figures  of  merit  could  occur;  for  example,  if  the  error  is 
ramping  from  positive  to  negative,  then  the  average  error  could  be 
zero  while  the  actual  error  varies  over  a  large  range  of  values. 
Evaluation  of  the  plots  of  mean  error  versus  time  provide  a  more 
accurate  and  complete  representation  of  tracker  performance. 

5.3  Performance  Plot  Format 

Performance  plots  are  generated  in  order  to  evaluate  fully  the 
tracking  algorithm’s  accuracy.  These  plots  are  of  x-  and  y-dynamics 
estimate  mean  errors  and  mean  error,  plus  and  minus  one  standard 
deviation  of  the  errors.  The  position  estimate  errors  in  pixels  are 
plotted  against  time  in  seconds.  Figure  5-1  is  the  y-position  error 
before  update  and  Figure  5-2  is  the  error  after  update.  The  maneuver 
takes  place  at  t  *  2.0  seconds;  the  plots  show  an  increase  in  error 
followed  by  a  recovery  to  approximately  zero-mean  tracking  error. 
Notice  that  the  maximum  mean  error  after  update  is  less  than  before 
update,  properly  showing  that  the  update  process  reduces  tracking 
error.  Figure  5-3  compares  the  actual  y-position  rms  error  with  the 
filter's  estimate  of  its  own  rms  error;  this  is  a  good  indication  of 
how  well  the  filter  is  tuned.  In  this  case,  the  maneuver  initiation 
causes  a  short  period  of  increased  error;  the  filter  responds 
appropriately  by  also  increasing  its  computed  rms  error.  This  plot 
format  is  used  to  tune  each  elemental  filter.  In  conjunction  with 
the  time  averaged  statistics,  the  mean  ±1  standard  deviation  error 
time  histories  provide  the  means  to  evaluate  the  tracking  algorithm's 
level  of  performance. 
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5.4  Parameter  Values 


The  actual  parameter  values  used  in  the  equations  developed  in 
earlier  chapters  can  significantly  affect  tracker  performance.  In 
all  cases,  the  values  used  for  this  research  are  selected  to  produce 
the  most  accurate  estimates  of  target  position.  Many  of  the  values 
are  the  result  of  previous  investigations  [8,13];  however,  some  are 
due  to  the  tuning  work  accomplished  during  this  effort.  Unless  so 
stated,  the  following  discussion  applies  equally  to  all  three 
filters . 

5.4.1  Truth  Model  Parameters.  The  initial  conditions  for  all 
of  the  trajectories  described  in  Section  3.2  are  identical  for  every 
simulation: 


Inertial  position: 


x  ■  5000  m 
y  *  500  m 
z 


20000  m 


Inertial  velocity: 


x  =  -1000  m/s 
y  =  0 

z  *  0 


Inertial  acceleration: 


x 

y 

z 


0 

0 

0 


The  maximum  intensity  of  each  hot  spot,  *max  of  Equation  (3-8)  is  20. 

The  variance  of  each  v^(t^)  in  Equation  (3-9),  accounting  for  both 

2 

FLIR  and  background  noises,  is  one  (intensity  unit)  .  This  leads  to 

a  signal-to-noise  ratio  of  20,  which  is  representative  of  many 

tracking  environments  [3].  Atmospheric  jitter  variance,  the  variance 

2 

of  xA  in  the  diagram  below  Equation  (3-1),  is  0.2  (pixels)  .  Based 
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on  a  glint  dispersion  parameter  of  2.0  (pixels)  ,  the  aspect  ratio, 
defined  in  Section  3.6,  is  one;  i.e.,  the  eigenvalues  of  in 
Equation  (3-8)  are  both  two.  These  parameters  applied  to  the  model 
developed  in  Chapter  III  specify  the  simulation  environment  used  to 
test  the  tracking  algorithms. 

5.4.2  Filter  Parameter  Values.  For  the  MMAFs,  each  element 
filter  was  tuned  to  produce  the  best  results  at  its  individual  design 
point  without  degrading  the  MMAF  performance.  In  the  two-element 
MMAF,  the  small  FOV  filter  required  a  high  dynamics  noise  variance  to 
prevent  MMAF  divergence  in  the  20g  scenarios  and  the  large  FOV  filter 
is  likewise  detuned  to  maintain  lock  on  a  constant  velocity  target. 
However,  the  three-element  MMAF  has  each  filter  tuned  closer  to  the 
optimum  values  for  the  type  of  trajectory  modeled  by  its  internal 
dynamics  model.  The  resulting  tuning  values  are  presented  in  Table 
IV- 1 . 

For  all  the  filters,  an  atmospheric  correlation  time  ( TAp)  of 

0.07  seconds  is  assumed.  The  filter  atmospheric  noise  variance  is 

2 

0.2  (pixels)  ,  like  the  truth  model. 

The  probability  lower  bound  is  set  at  0.001.  The  initial 
probability  of  filter  one  is  0.99  for  the  two  element  MMAF  and  0.98 
for  the  three-element  MMAF.  This  high  value  indicates  confidence 
that  the  target  is  essentially  at  constant  velocity  and  not 
performing  evasive  maneuvers. 


Although  earlier  chapters  discuss  every  aspect  of  the  tracking 
algorithms  being  investigated  and  the  target  test  trajectories,  a 
brief  description  of  each  will  enhance  understanding  of  the 
performance  analysis  contained  in  the  next  chapter. 

5.5.1  Multiple  Model  Adaptive  Filters.  Three  different 
tracking  algorithms  are  being  evaluated,  the  two-element  baseline 
MMAF,  the  three-element  Bayesian  MMAF,  and  the  three-element  MAP 
MMAF.  The  baseline  MMAF  is  taken  directly  from  Suizu's  work  [13] 
based  on  the  Millner  algorithm  [8]  for  the  elemental  filters.  In 
addition  to  providing  a  basis  for  evaluating  the  three-element 
filters,  the  baseline  MMAF  is  exercised  for  eight  seconds  (rather 
than  five  seconds)  and  operated  at  100  Hz.  These  variations 
investigate  crossover  effects  and  the  performance  benefit  of 
increased  sample  rate.  The  latter  two  MMAFs  contain  the  same 
elemental  filters;  however,  they  create  the  MMAF  estimate 
differently.  The  Bayesian  MMAF  produces  the  weighted  average  of  the 
elemental  filter's  estimates,  using  the  computed  conditional 
probabilities  as  developed  in  Chapter  II.  Also  using  the  conditional 
probabilities,  the  MAP  MMAF  chooses  the  elemental  estimate  with  the 
highest  probability  of  being  correct.  The  performance  of  these  three 
filters  is  evaluated  over  a  number  of  test  cases. 

5.5.2  Test  Cases.  Each  of  the  MMAFs  are  exercised  against  the 
same  combination  of  trajectories  at  different  acceleration  levels,  as 
shown  in  Table  5.1  below. 


Table  5.1.  Test  Cases 


BASELINE 

BAYESIAN 

MAP 

Traj  1  -  straight  line 

X 

X 

X 

Traj  2  -  pull-up 

g  level: 

2 

X 

X 

X 

10 

X 

X 

X 

20 

X 

X 

X 

Traj  3  -  pull-up/ 

constant  velocity 

climb 

g  level: 

2 

X 

X 

X 

10 

X 

X 

X 

20 

X 

X 

X 

Trajectory  four  is  not  used  because  it  is  designed  to  investigate 
target  image  estimation  rather  than  tracking  accuracy.  The  target 
remains  in  approximately  the  same  position  in  the  FOV,  changing  only 
its  range  and  orientation  and  does  not  stress  the  algorithm's  ability 
to  estimate  target  position.  Additionally,  each  elemental  filter  is 
tested  against  trajectory  two  at  the  g  level  for  which  it  is  tuned. 
The  elemental  filter's  performance  with  optimum  tuning  should 
represent  the  "best"  tracking  performance  in  that  test  case;  because 
of  purposeful  detuning  of  various  degrees,  the  elemental  filter's 
performance  may  not  be  the  "best".  However,  the  comparison  between 
the  elemental  filter  and  the  MMAF  gives  some  idea  of  the  effect  the 
multiple  model  approach  has  on  tracking  performance. 

A  test  case  designation  code  is  useful  to  identify  the  different 
combinations  of  filters,  trajectories,  and  g  levels.  The  following 
code  is  used  to  identify  the  test  cases  in  tables  and  plots. 


65 


T2 


G20 


F2- 


BAY 


1 


2  3 


4 


1.  Trajectory  number 

2.  Maneuver  g  level  (2,10,20) 


3. 


Filter  Type 
MM  -  multiple  model 


FI 

F2  V  element 
F3 


filters 


4.  Filter  variation 

BAY  -  Bayesian  MMAF 

MAP  -  MAP  MMAF 

BASE  -  Baseline  MMAF 

XT  -  Extended  time,  baseline  MMAF 

100HZ  -  100Hz  sample  rate,  baseline  MMAF 


Under  this  scheme,  the  example  code  (T2G20F2-BAY)  indicates  the 
second  elemental  filter  in  the  Bayesian  MMAF  operating  in  a  20g, 
trajectory  two. 

5.6  Summary 

The  evaluation  methods  presented  in  this  chapter  are  used  to 
describe  the  results  obtained  from  operating  the  three  MMAFs  under 
various  conditions.  These  conditions  are  described  by  the  selected 
target  trajectories  and  parameter  values  for  the  truth  model  and  the 
filters.  Based  on  the  statistics  and  plot  formats  discussed  in  this 
chapter,  the  next  chapter  investigates  the  performance  of  the 
different  MMAF  formulations. 


66 


1  VI.  PERFORMANCE  ANALYSIS 

6. 1  Introduction 

•  The  performance  analysis  of  the  three  MMAFs  is  discussed  in  this 

chapter.  Each  filter  is  evaluated  separately,  based  on  both  its 
temporal  averages  and  error  statistics  times  histories.  The  baseline 
1  MMAF  is  presented  first,  along  with  the  extended  simulation  (8  sec) 

and  the  100  Hz  sample  rate  version  of  this  filter.  The  addition  of  a 
third  elemental  filter  is  investigated  next  by  analyzing  the  Bayesian 
MMAF.  Finally,  the  MAP  MMAF  results  are  evaluated.  The  concluding 
section  compares  and  contrasts  the  results  of  all  three  MMAFs. 

j  %  *  6.2  The  Baseline  MMAF 

As  discussed  in  Section  4.4,  the  narrow-FOV  elemental  filter 

requires  a  high  level  of  dynamics  driving  noise  to  prevent  total  loss 

|  of  target  at  20  g's  by  the  baseline  MMAF.  For  all  trajectories  at 

low  accelerations  (1-2  g's),  the  conditional  probabilities  maintain 

a  heavy  weight  (0.99+)  on  filter  1.  This  indicates  that  the  actual 

T 

residuals  remain  well  matched  to  the  internally  computed  (H  £  II  +  R) 
of  that  elemental  filter.  The  high  level  of  acceleration  noise  on 
the  dynamic  states  allows  this  model  to  track  a  range  of  lower 
acceleration  values  rather  than  just  zero  acceleration. 

In  trajectory  two,  significant  variations  of  the  conditional 
probabilities  occur  at  10  and  20  g's.  At  10  g's,  the  shift  of 


weighting  to  filter  2  starts  within  0.3  seconds  after  the  maneuver  is 


initiated,  and  after  two  thirds  of  a  second  the  weighting  is  heavily 
on  filter  2.  The  MMAF  responds  slightly  faster  at  20  g's.  No 
reacquisitions  occur  at  10  g’s;  however,  the  20  g  maneuver  averages 
two  reacquisitions  per  run  for  the  narrow  FOV  filter.  This  is 
indicative  of  the  narrow  FOV  filter  having  severe  difficulty  in 
maintaining  lock  when  a  20  g  step  acceleration  is  introduced.  Use  of 
a  more  realistic  acceleration  model  might  remove  the  reacquisition 
requirement  altogether.  The  conditional  probabilities,  following  the 
maneuver,  tend  to  oscillate  in  the  0.4  to  0.5  range  and  sometimes 
reach  as  high  as  0.7.  This  is  due  to  the  purposeful  detuning  of  the 
elemental  filters  to  avoid  total  loss  of  the  target  at  20  g's; 
neither  elemental  filter  produces  residuals  that  are  significantly 
large  enough  to  push  the  weighting  one  way  or  the  other.  Up  until 
the  3.5  second  point  in  the  simulation,  trajectory  three  behavior  is 
the  same;  however,  at  this  point  the  maneuver  change  in  this 
trajectory  causes  the  wide  FOV  to  be  weighted  much  more  heavily  than 
the  narrow  FOV.  Once  again,  as  the  target  trajectory  remains 
consistent  after  this  step  change  in  acceleration  occurs,  the  MMAF 
begins  to  increase  the  weighting  on  filter  i. 

The  plots  of  mean  error  ±  one  standard  deviation  in  the 
x-position,  for  trajectory  two  (Appendix  A),  do  not  show  significant 
changes  due  to  the  maneuver,  because  the  actual  trajectory  change  in 
the  x-direction  is  relatively  small.  The  x-position  mean  error 
begins  a  positive  ramp  after  the  maneuver  at  2.0  seconds;  the  maximum 
rate  is  0.5  pixels/second  during  the  20  g  maneuver.  This  ramping  can 
possibly  be  attributed  to  the  simplifying  assumptions  incorporated  in 


the  models.  First,  each  elemental  filter  assumes  a  non-rotating 
inertial  reference  frame  when,  in  fact,  the  tracker-based  coordinate 
frame  does  rotate.  This  model  simplification  can  be  mitigated  by  the 
addition  of  pseudo-noise  to  the  position  and  velocity  state 
estimates.  Second,  the  true  position  is  the  location  of  the  target 
center  of  mass,  and  the  filter  is  tracking  the  target  radiation 
centroid;  usually  the  two  are  not  coincident.  A  definitive 
explanation  of  the  ramping  has  not  yet  been  established,  and  this 
warrants  further  investigation. 

Meanwhile,  the  y-position  has  a  sharp  increase  in  error 
following  the  maneuver,  but  quickly  recovers  (0.625  seconds)  to  zero- 
mean  tracking  error.  The  one  sigma  errors  are  slightly  higher  after 
the  maneuver  due  to  the  increased  emphasis  on  filter  2.  The  time 
averaged  statistics  in  Table  6-1,  show  this  and  also,  the  tracking 
error  is  reduced  following  updates.  Additionally,  Table  6-1 
indicates  the  effects  of  multiple  model  weighting  on  filter 
performance.  Comparing  the  MMAF  and  filter  2  in  a  20  g  trajectory 
two,  the  single  filter  has  consistently  better  performance.  The  same 
type  of  comparison  for  trajectory  one  indicates  that  filter  1  does 
not  perform  noticeably  different  from  the  MMAF,  probably  due  to  the 
consistently  high  weighting  on  filter  1  in  the  MMAF.  Although  it  is 
expected  that  the  tracking  error  will  increase  slightly  for  higher  g 
levels,  one  notes  that  the  MMAF  actually  performs  worse  in  the  well- 
modeled  20  g  maneuver  than  in  the  10  g  maneuver  that  is  not 
explicitly  modeled  by  one  of  the  elemental  filters.  This  is  further 
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evidence  of  the  purposeful  detuning  of  the  elemental  filters  which 
leads  to  the  oscillation  in  weighting  of  the  two  elemental  filters 
discussed  earlier  in  this  section. 

Table  6-1.  Time  Averaged  Statistics  (mean  ±  1  o ) 
for  the  Baseline  MMAF  (3.5  -  5.0  sec) 


CASE 

x(-) 

(pixels) 

x(+) 

(pixels) 

y(-) 

(pixels) 

(y+) 

(pixels) 

TIMM-BASE 
(T  IF  1-BASE) 

232 

.216 

+  .424 
±.426 

-.219 

.203 

±.359 

±.377 

-.000 

-.009 

±.394 

±.394 

-.003 

-.011 

±.345 

±.345 

T2G2MM-BASE 

.214 

±.424 

.198 

±.369 

-.324 

±.404 

-.008 

±.355 

T2G10MM-BASE 

.465 

±.431 

.389 

±.374 

-.120 

±.451 

.008 

±.402 

T2G20MM-BASE 

(T2G20F2-BASE) 

1.29 

.216 

±.444 

±.420 

1.04 

.203 

±.388 

±.365 

-.179 

-.009 

±.467 

±.394 

.044 

-.011 

±.422 

±.345 

T3G2MM-BASE 

.127 

±.481 

.158 

±.406 

.181 

±.405 

.131 

±.358 

T3G10MM-BASE 

.119 

±.516 

.157 

±.429 

.553 

±.506 

.420 

±.441 

T3G20MM-BASE 

.305 

±.501 

.376 

±.415 

1.07 

±.514 

.816 

±.441 

The  error  statistics  time  history  plots  for  x-position  on 
trajectory  three  (Appendix  A)  are  significantly  different  from 
trajectory  two  plots.  At  approximately  1.5  second  into  the 
simulation,  the  mean  error  experiences  a  step  change,  causing  a  bias 
until  the  maneuver  at  3.5  seconds.  The  bias  value  is  approximately 
1.5  pixels  for  all  acceleration  levels;  rms  error  levels  remain  about 


the  same.  Returning  to  approximately  zero-mean  error  at  3.5  seconds, 
the  x-position  plot  demonstrates  the  characteristic  ramping  mean 
error  as  seen  in  trajectory  two.  There  are  no  known  differences 
between  trajectories  two  and  three  in  the  first  3.5  seconds  of  the 
simulation  to  account  for  this  unexpected  behavior;  it  is  supposed 
that  the  trajectory  generation  in  the  truth  model  may  have  an 
undiscovered  variation. 

Tracking  performance  in  the  y-direction  is  considerably  better 
than  the  x-direction.  This  result  is  highly  encouraging  because 
tracking  in  the  y-direction  is  significantly  more  difficult  than  in 
the  x-direction;  target  maneuvers  show  the  greatest  change  in 
y-direction  motion.  Furthermore,  errors  in  position  estimation  along 
the  length  of  the  target  are  not  as  serious  as  errors  orthogonal  to 
the  target  longitudinal  axis.  The  maneuver  changes  result  in  a 
maximum  mean  excursion  of  -5.25  pixels  as  2.0  seconds,  followed  by  a 
maximum  of  2.0  pixels  after  the  maneuver  at  3.5  seconds.  The  filter 
recovers  in  approximately  0.75  seconds  following  the  first  maneuver; 
the  response  to  the  second  change  is  somewhat  slower,  possibly  due  to 
the  higher  weighting  of  the  wide  FOV  elemental  filter.  The  time 
averages  in  Table  6-1  are  artificially  high  because  the  averaging  is 
performed  from  the  3.5  second  maneuver  point  until  simulation  end; 
the  plots  also  show  a  slow  return  to  zero-mean  tracking  error. 

6.2.1  Extended  Time  Baseline  MMAF.  The  baseline  MMAF,  as 
implemented  by  Suizu  [13],  demonstrated  a  tendency  to  increase  mean 
tracking  error  at  approximately  5.0  seconds.  This  time  represents 
the  minimum  range/maximum  cross  rate  condition  for  the  simulation;  to 
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investigate  tracking  behavior  beyond  this  condition,  the  simulation 
time  was  extended  to  8.0  seconds.  In  most  of  the  test  cases,  the 
extended  runs  demonstrated  the  same  behavior  as  the  baseline  MMAF 
with  a  few  notable  exceptions. 

The  conditional  probability  behavior  is  virtually  identical  to 
the  baseline  MMAF  for  maneuvers  of  1  -  10  g's.  However,  despite 
similar  trends  in  system  identification  (e.g.,  appropriate  heavy  p^ 
values  on  the  best  matched  filter) ,  the  extended  runs  show  the  MMAF 
losing  lock  in  20  g  maneuvers.  Unlike  the  same  MMAF  in  a  5.0  second 
run,  loss  of  lock  occurs  within  0.5  seconds  of  the  maneuver  at  2.0 
seconds.  This  failure  to  maintain  lock  is  attributed  to  a  different 
noise  time  history  caused  by  additional  calls  to  the  random  number 
generator  used  to  create  the  additive  noise  values  for  the  additional 
frames  of  data  in  the  extended  runs.  This  loss  of  lock  indicates 
that  the  tuning  values  are  very  marginal  in  the  baseline  MMAF. 

The  x-position  mean  ±  one  standard  deviation  error  plots 
(Appendix  D)  for  trajectory  three  differ  significantly  from  the 
baseline  MMAF  plots.  In  this  case,  the  unexplained  step  bias  at  1.5 
seconds  is  not  observed.  When  the  filter  maintains  lock,  the 
y-position  plots  demonstrate  recovery  at  the  same  rate  as  the 
baseline  MMAF.  These  trends  are  further  illustrated  by  the  time 
average  statistics  in  Table  6-2.  As  seen  in  the  table,  these 
statistics  show  slightly  better  mean  tracking  errors  than  the 
baseline,  however,  the  values  are  averaged  later  in  the  simulation 
after  the  t  .Liter  has  more  time  to  recover  from  maneuver  transients. 
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Table  6-2.  Time  Averaged  Statistics  (mean  i  1  a) 


for  the 

Extended  Runs 

(6-8  sec) 

x(-) 

x(+) 

y(-) 

(y+) 

CASE 

(pixels) 

(pixels) 

(pixels) 

(pixels) 

TIMM- XT 

.405  ±.436 

.373  +.385 

.005  +.426 

.001  ±.378 

T2G2MM-XT 

.410  +.433 

.375  +.385 

-.020  +.427 

-.003  ±.377 

T2G10MM-XT 

.800  +.459 

.685  +.412 

.042  +.424 

.115  ±.373 

T3G2MM-XT 

.449  +.432 

.418  +.382 

.064  +.418 

.060  ±.369 

T3G10MM-XT 

.397  +.429 

.367  +.377 

.143  +.419 

.132  +.370 

6.2.2 

Baseline  MMAF  - 

100  H.z  Sample 

Rate.  Optical 

implementa- 

tion  of  the  tracking  algorithm,  as  proposed  by  Roemer  [10]  and  Rogers 
[11],  allows  data  to  be  processed  at  the  increased  sample  rate  of  100 
Hz.  All  other  aspects  remaining  constant,  additional  updates  should 
greatly  improve  tracking  accuracy.  In  actuality,  increasing  the 
baseline  30  Hz  sample  rate  to  100  Hz  did  not  improve  tracking 
performance.  Under  this  condition,  the  MMAF  diverged  in  every  case 
and  failed  to  maintain  lock  at  higher  (10  -  20  g's)  acceleration 
levels.*  This  type  of  behavior  may  indicate  problems  in  the 
measurement  model.  The  increase  in  sample  rate  reduces  the  smoothing 
effect  of  the  propagation  portion  of  the  algorithm,  allowing 
mismodeled  measurement  data  to  have  a  more  severe  effect  on  the  state 
estimate.  A  poor  measurement  model  could  lead  the  filter  into 
erroneous  evaluations  of  the  incoming  data  causing  poor  tracking 
performance.  Due  to  limitations  on  time  and  computer  resources,  this 

*  Subsequent  investigation  located  a  coding  error.  See  Appendix  E 
for  updated  results. 
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difficulty  was  not  fully  investigated;  however,  the  plots  which  did 
not  demonstrate  loss  of  lock  are  included  in  Appendix  E  for 
completeness . 

6.3  Bayesian  MMAF 

The  Bayesian  MMAF  consists  of  the  three  elemental  filters 
described  in  Section  4.4.  This  filter  is  tested  against  the  same 
trajectories  as  the  baseline  MMAF.  In  all  of  the  trajectories  at  low 
accelerations  (1  -  2  g's),  the  rate  of  change  of  the  conditional 
probabilities,  following  a  maneuver,  are  very  similar  to  the  baseline 
MMAF;  at  the  higher  accelerations,  the  trends  differ  markedly.  At  10 
g's  in  trajectory  two,  following  the  maneuver,  filter  2  is  weighted 
more  heavily  than  the  other  two  filters,  but  within  1  second  filter  1 
is  once  again  strongly  weighted.  In  the  20  g  maneuver,  the  weighting 
shifts  to  filter  3  for  a  few  frames  and  then  to  filter  2,  followed 
immediately  by  the  narrow  FOV  filters  (1  and  3)  losing  lock.  After 
reacquisition,  the  conditional  probability  of  filter  2  is  very  high 
(0.9+)  and  remains  so  until  1.3  seconds  when  filter  1  rapidly 
dominates  the  weighted  average.  In  each  case,  less  of  the 
oscillations  in  probabilities  previously  observed  in  the  baseline 
MMAF  occurred.  As  observed  in  the  baseline  MMAF,  when  the  filter 
obtains  good  acceleration  estimates,  the  weighting  returns  to  filter 
1.  Trajectory  three  also  has  the  same  conditional  probability  trends 
as  observed  in  the  baseline  cases.  However,  following  the  maneuver 
change  at  3.5  seconds,  the  filter  is  much  slower  in  returning 
emphasis  to  filter  1;  filter  l’s  probability,  in  a  10  g  acceleration, 


only  reaches  0.68  as  compared  to  0.9+  in  the  baseline  case,  and  the 
20  case  fails  to  maintain  lock.  Loss  of  lock  occurs  very  quickly 
following  the  maneuver  at  2.0  seconds;  this  is  an  additional 
unidentified  peculiarity  in  trajectory  three  (in  that  it  should  yield 
the  same  characteristics  as  trajectory  two  until  the  3.5  second 
point) . 

Tracking  performance  in  the  x-direction  remains  fairly 
consistent,  including  the  unexplained  bias  observed  only  in 
trajectory  three  (Appendix  B) .  The  time  averaged  statistics  in  Table 
6-3  show  that  the  mean  error  tends  to  increase  at  higher 
accelerations,  as  does  the  error  standard  deviation.  In  most  cases, 
the  elemental  filters,  operating  at  the  conditions  for  which  they  are 
specifically  tuned,  demonstrate  lower  mean  errors  and  higher  standard 
deviations  than  the  MMAF  in  the  same  test  case.  This  results  from 
the  MMAF  mean  value  not  being  corrupted  by  the  inclusion  of  means 
from  mismatched  filters;  similarly,  the  covariance  values  are  higher 
because  the  values  in  filters  2  and  3  are  not  reduced  by  being 
combined  with  the  lower  value  used  in  filter  1. 

Plots  of  error  statistics  time  history  demonstrate  very  good 
performance  in  the  y-direction  (Appendix  B) .  In  every  case  the  MMAF 
recovers  within  0.75  seconds  following  a  maneuver;  this  recovery  time 
is  the  same  for  the  elemental  filters.  The  elemental  filter  that  is 
best  tuned  to  the  given  maneuver  (filter  2  or  3) ,  alone,  demonstrates 
lower  maximum  excursions  following  a  maneuver  than  the  MMAF. 
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Table  6-3.  Time  Averaged  Statistics  (mean  £  la) 

for  the  Bayesian  MMAF  (3.5  -  5.0  seconds) 


I 


CASE 

x(-) 

(pixels) 

x(+) 

(pixels) 

y(-) 

(pixels) 

(y+) 

(pixels) 

TIMM-BAY 

.228 

±.  422 

.215 

£.367 

.026  £.419 

.025  £.372 

(T IF  1-BAY) 

.216 

£.420 

.203 

£.365 

.009  £.394 

-.012  £.345 

T2G2MM-BAY 

.220 

£.429 

.205 

£.374 

-.012  £.427 

.013  £.378 

T2G10MM-BAY 

.450 

£.456 

.376 

£.400 

-.239  £.458 

-.101  £.406 

(T2G10F3-BAY) 

.388 

£.482 

.321 

£.411 

-.420  £.429 

-.245  £.368 

T2G20MM-BAY 

1.13 

£.465 

1.10 

£.411 

-.271  £.481 

-.101  £.44? 

(TG20F2-BAY) 

1.02 

£.482 

.793 

£.414 

.062  £.450 

.331  £.378 

T3G2MM-BAY 

.132 

£.507 

.170 

£.422 

.130  £.434 

.090  £.376 

T3G10MM-BAY 

.088 

£.511 

.129 

£.427 

.367  1.440 

.226  £.340 

T3G20MM-BAY 

lost 

lock  on  target 

However,  these  two  elemental  filters  demonstrate  non-zero  mean  trends 
following  recovery;  filter  2  has  a  positively  ramping  bias  and  filter 
3  never  quite  returns  to  zero-mean  tracking  error:  it  maintains  a 
small  negative  bias  until  the  end  of  the  simulation.  The  MMAF  does 
recover  to  zero-mean  tracking;  although  in  trajectory  three  at  20 
g's,  it  recovers  more  slowly  than  at  10  g's,  probably  due  to  the 
longer  emphasis  on  the  wide  FOV  filter.  Generally,  the  Bayesian  MMAF 
outperforms  the  baseline  MMAF;  this  is  discussed  in  more  detail  in 


Section  6.5 


6.4  MAP  MMAF 


This  variation  of  the  MMAF  uses  the  output  of  the  single  "best" 
elemental  filter  to  form  the  MMAF  estimate.  The  selection  is  made 
based  on  the  calculated  conditional  probabilities.  In  the  lower 
acceleration  cases,  the  MAP  algorithm  selects  filter  1  exclusively 
and  the  actual  calculated  probabilities  for  filter  1  are  higher  than 
those  observed  for  filter  1  in  the  Bayesian  MMAF.  Possibly, 
reacquisition,  as  discussed  in  Section  4.4,  allows  the  diverging 
filters  to  be  inappropriately  weighted  for  a  few  frames,  thus  leading 
to  different  conditional  probability  weights  subsequently.  Following 
the  10  g  maneuver  in  trajectory  two,  filter  2  is  selected  within  0.3 
seconds;  for  the  next  second,  the  selection  rotates  through  3,  2,  1, 
and  3  in  order.  After  4.16  seconds,  filter  1  is  used  for  the 
remainder  of  the  simulation.  Similar  behavior,  about  2  sample 
periods  sooner,  is  observed  at  20  g's;  the  narrow  FOV  filters  undergo 
reacquisition  at  about  2.5  seconds.  Trajectory  three,  at  both  10  and 
20  g's,  demonstrates  the  same  trends  until  3.5  seconds  when  the 
second  maneuver  change  occurs.  Following  the  return  to  constant 
velocity  flight,  the  MAP  MMAF  begins  to  shift  between  filters  2  and 
3.  Filter  1  is  selected  again  by  4.67  seconds,  at  10  g's;  however, 
at  20  g's,  filters  2  and  3  continue  to  alternate  until  the  end  of  the 
simulation.  Thus,  as  observed  previously,  changes  in  acceleration 
trigger  a  shift  to  filters  2  and  3,  while  filter  1  is  selected  when  a 
good  estimate  of  constant  acceleration  i>  obtained  [1], 


Tracking  performance  in  the  x-direction  is  consistent  with  that 
in  the  other  two  MMAFs  (Appendix  C) .  Generally,  as  shown  in  Table 
6-4,  the  MAP  MMAF  has  lower  mean  errors  and  higher  error  standard 
deviations  than  the  Bayesian  MMAF.  The  ramp  bias  observed  in  both 
the  baseline  and  the  Bayesian  MMAFs  still  occurs  with  approximately 
the  same  magnitude.  Additionally,  the  unresolved  step  bias  observed 
in  trajectory  three  still  remains;  however,  the  bias  magnitude  is 
slightly  higher,  1.6  and  1.9  pixels  for  10  and  20  g's,  respectively, 
as  compared  to  1.5  pixels  in  the  first  two  MMAFs.  The  exclusive  use 
of  filter  1  data  seems  to  make  the  output  more  sensitive  to  this 
phenomena. 

Table  6-4.  Time  Averaged  Statistics  (mean  ±  la) 
for  the  MAP  MMAF  (3.5  -  5  sec) 


CASE 

x(-) 

(pixels) 

x(+) 

(pixels) 

y(-) 

(pixels) 

(y+) 

(pixels) 

T1MM-MAP 

.216 

±.420 

.203 

±.366 

.005 

±.416 

.003 

±.350 

T2G2MM-MAP 

.216 

±.428 

.198 

±.375 

-.014 

±.432 

.010 

±.387 

T2G10MM-MAP 

T2G10F3-MAP* 

.453 

.400 

±.446 

±.468 

.377 

.333 

±.391 

±.468 

-.260 

-.311 

±.528 

±.608 

-.129 

-.126 

±.481 

±.547 

T2G20MM-MAP 

1.26 

±.461 

1.00 

±.412 

-.344 

±.478 

-.129 

±.454 

T3G2MM-MAP 

.109 

±.520 

.149 

±.431 

.136 

±.443 

.095 

±.381 

T3G10MM-MAP 

.075 

±.512 

.116 

±.427 

.359 

±.454 

.219 

±.383 

T3G20MM-MAP 

.159 

±.501 

.247 

±.416 

.991 

±.525 

.697 

±.446 

*  The  dynamics  variances  on  target  acceleration  for  the  n 
filters  1,  2,  and  3  are  300,  1900,  and  2000  pixels  /sec 


As  might  be  expected,  the  y-position  error  statistics  plots 
(Appendix  C)  demonstrate  the  effects  of  MAP  estimation  more  clearly. 
Low  acceleration  (0  -  2  g's)  maneuvers  are  tracked  entirely  by  filter 
1.  The  maximum  excursion  of  mean  error  observed  is  0.5  pixels  and 
the  error  standard  deviation  remains  under  0.5  pixels.  In  this  case, 
accuracy  has  been  sacrificed  by  the  purposeful  detuning  of  filter  1, 
required  for  the  MMAF  to  maintain  lock  in  20  g  maneuvers.  Trajectory 
two  shows  a  slight  tendency  toward  a  negative  bias  following  recovery 
from  the  maneuver.  At  10  and  20  g's,  the  maximum  excursions  are  -1.6 
and  -4.0  pixels  and  the  recover  time  is  0.68  and  0.75  seconds, 
respectively.  The  increased  standard  deviation  due  to  the  selection 
of  filter  2  or  3,  following  the  maneuver,  is  very  noticeable  in  the 
10  g  case  and  less  so,  due  to  scaling,  in  the  20  g  case.  In  an 
effort  to  evaluate  the  effect  of  detuning,  the  elemental  filters  are 
tuned  more  tightly,  using  the  values  noted  in  Table  6-4  rather  than 
those  established  in  Section  4.4.  The  MAP  MMAF  is  tested  at  10  g's 
(loss  of  lock  by  the  MMAF  at  20  g's  is  the  only  reason  for  detuning 
the  elemental  filters);  the  MMAF  with  this  tuning  has  a  slightly 
smaller  maximum  excursion  (-1.5  pixels)  and  quicker  recovery  (0.625 
seconds) ,  but  it  shows  a  greater  tendency  towards  a  negative  bias 
following  recovery  (Appendix  C) .  The  time  averaged  statistics,  in 
Table  6-4,  indicate  that  most  of  the  improvement  is  in  the 
x-direction  tracking  rather  than  the  y-direction  which,  as  discussed 
earlier,  is  not  the  difficult  direction  to  track.  In  this  light,  the 
current  tuning  values  appear  more  favorable  than  might  be 
anticipated.  Trajectory  three  performance  also  demonstrates  these 


trends.  The  slightly  higher  standard  deviations  following  maneuvers 
are  observable  but  not  obvious  in  the  plots.  The  MAP  estimates  are 
not  significantly  lower  in  this  case  than  the  other  tow  filters;  this 
appears  to  be  a  consequence  of  the  cycling  between  filters  which 
leads  to  an  effective  averaging  of  the  elemental  filter's  outputs. 

6.5  Comparison  of  the  Three  MMAFs 

A  comparison  of  the  three  different  MMAFs  (baseline,  Bayesian, 
and  MAP)  allows  one  to  evaluate  the  effects  of  the  third  filter  in 
the  bank  and  the  MAP  estimation  approach  versus  the  Bayesian 
approach.  Several  aspects  of  each  filter  are  used  to  judge  relative 
performance.  The  maximum  excursion  of  the  mean  error  following  a 
maneuver  indicates  how  well  the  MMAFs  handle  changing  target  motion. 
Additionally,  the  time  to  recover  to  zero  mean  tracking  error  after 
the  maneuver  is  another  measure  of  performance  in  this  category. 
Observation  of  the  time  averaged  statistics  yields  an  idea  of  each 
MMAF's  accuracy  in  steady  state. 

Clearly,  as  seen  in  Table  6-5,  the  addition  of  the  third  filter 
to  the  bank  reduces  the  maximum  mean  error  following  the  maneuvers  at 
2.0  and  3.5  seconds.  The  differences  between  the  Bayesian  and  MAP 
MMAFs  do  not  indicate  a  significant  difference  between  the  two 
approaches;  however,  the  MAP  MMAF  may  handle  20  g  maneuvers  better. 
The  evidence  is  inconclusive  due  to  the  difficulties  with  trajectory 
three  discussed  earlier  in  this  chapter. 
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Table  6-5.  Maximum  Mean  Error  Excursion  Following 
Maneuvers  at  2.0  and  3.5  Seconds 


Trajectory  Two 

(2.0  sec) 

Trajectory  Three 

(2. 0/3/5  sec) 

10  g's 

20  g's 

10  g's 

20  g's 

(pixels) 

(pixels) 

(pixels) 

(pixels) 

Baseline 

-2.38 

-6.5 

-2.38/1.25 

-6. 5/2.0 

Bayesian 

-1.4 

-4.13 

-1.51/1.18 

lost  lock 

MAP 

-1.6 

-4.00 

-1.51/1.11 

-5.38/2.5 

Response  time  after  a  maneuver  (Table  6-6)  doesn't  vary 
significantly  between  the  three  MMAFs.  Immediately  following  a 
maneuver,  when  the  conditional  probabilities  shift,  it  takes  about 
the  same  number  of  sample  periods  to  cause  a  change  of  filter  in  the 
MAP  MMAF  as  to  achieve  a  0.9+  weighting  of  the  same  filter  in  the 
Bayesian  MMAF. 


Table  6-6.  Recovery  Time  After  Maneuvers  at 
2.0  and  3.5  Seconds 


Trajectory  Two 

(2.0  sec) 

Trajectory  Three 

(2. 0/3/5  sec) 

10  g's 

20  g's 

10  g's 

20  g's 

(sec) 

(sec) 

(sec) 

(sec) 

Baseline 

.625 

.75 

.688/  * 

.75/  * 

Bayesian 

.75 

.75 

.688/. 75 

lost  lock 

MAP 

.688 

.75 

.688/. 688 

.813/. 938 

*  Never  reached  zero-mean  tracking  error  in  the  5.0  sec 


The  time  averaged  statistics  shown  in  Table  6-7  do  not  indicates 
any  overwhelming  trends.  There  is  no  strong  evidence  that  the 
addition  of  the  third  filter  improves  steady  state  tracking  error  in 
trajectory  two;  however,  trajectory  three,  y-direction,  demonstrates 
reduced  mean  errors  except  when  the  Bayesian  MMAF  loses  lock.  The 
MAP  MMAF  does  demonstrate  generally  lower  means  and  higher  standard 
deviations  than  the  Bayesian  MMAF,  but  not  significantly.  Notice 
that  the  standard  deviations  remain  consistently  between  about  0.35 
and  0.45;  this  is  very  important  for  insuring  that  the  laser  energy 
"paints"  a  very  small  area  of  the  target. 


Table  6-7.  Time  Averaged  Statistics  (mean  ±  la) 
After  Update 


Trajectory  Two 

Trajectory  Three 

x(+) 

y(+) 

x(+) 

y(+) 

(pixels) 

(pixels) 

(pixels) 

(pixels) 

2  g's 

Baseline 

.198  ±.369 

.008  ±.355 

.158 

± .  406 

.131  ±.358 

Bayesian 

.205  ±.374 

.013  ±.378 

.170 

±.422 

.090  ±.376 

MAP 

.198  ±.375 

.010  ±.387 

.149 

±.431 

.095  ±.381 

10  g’s 

Baseline 

.389  ±.374 

.008  ±.422 

.157 

±.429 

.420  ±.441 

Bayesian 

-.376  ±.400 

-.101  ±.406 

.129 

±.427 

-.226  ±.340 

MAP 

.377  ±.391 

-.129  ±.481 

.116 

±.427 

.219  ±.383 

20  g's 

Baseline 

1.04  ±.376 

.044  ±.422 

.376 

±.415 

.816  ±.441 

Bayesian 

-1.13  ±.456 

.065  ±.449 

lost 

lock 

MAP 

1.00  ±.412 

-.129  ±.457 

.247 

±.416 

.697  ±.446 

The  ability  to  track  a  highly  dynamic  target  is  enhanced  by  the 
addition  of  a  third  elemental  filter  to  the  multiple  model  adaptive 
filter  bank.  Lower  peak  mean  errors  increase  the  likelihood  of 
achieving  a  kill  with  a  laser  weapon  system.  The  additional 
computational  loading  required  by  a  third  element  becomes  almost 
insignificant  if  parallel  processing  techniques  are  employed. 

The  baseline  MMAF  does  not  encounter  any  difficulty  with  the 
minimum  range/maximum  crossing  rate  condition;  however,  it 
demonstrates  divergent  characteristics  at  a  100  Hz  sample  rate.  This 
poor  performance  at  an  increased  sample  rate  indicates  a  need  to 
question  the  basic  measurement  model  employed  in  each  elemental 
filter. 

Initially,  it  was  expected  that  the  differences  between  Bayesian 
and  MAP  estimation  would  be  more  noticeable.  The  anticipated  faster 
response  of  the  MAP  MMAF  to  maneuvers  was  not  demonstrated.  Average 
mean  tracking  errors  and  the  maximum  mean  excursions  were  very 
similar  in  both  approaches.  The  inability  of  the  algorithms  to 
maintain  consistent  system  identification,  e.g.,  maintain  high 
conditional  probabilities  on  the  "best"  matched  model,  could  explain 
these  ambiguous  results.  Additionally,  the  reacquisition  function 
may  have  prevented  more  accurate  and  consistent  system 
identification;  however,  because  it  only  appeared  during  20  g 
maneuvers,  it  is  a  relatively  minor  influence.  Overall,  the  study 


shows  that  an  additional  element  in  the  MMAF  bank  improves 
performance  and  that  both  Bayesian  and  MAP  estimation  techniques 
support  accurate  tracking  of  highly  dynamic  airborne  targets. 


7.2  Recommendations 

Further  investigation  of  the  multiple  model  adaptive  filter  to 
resolve  the  difficulties  encountered  and  expand  current  knowledge  is 
recommended.  Future  efforts  should  address  the  following  areas: 


Development  of  more  realistic  target  trajectories. 
Specifically,  the  step  changes  in  accelerations  should  be 
replaced  with  models  that  more  closely  approximate  typical 
target  behavior.  This  may  allow  the  MMAF  to  maintain  lock  with 
tighter  tuning  of  the  elemental  filters. 

Investigate  the  step  bias  in  x-position  observed  only  in 
trajectory  three 

Test  the  tracking  algorithms  against  targets  with  varying  shape 
functions  caused  by  changes  in  target  orientation  and  variations 
in  size,  shape,  and  number  of  hot  spots 

Indepth  investigation  of  the  100  Hz  sample  rate,  concentrating 
on  possible  inaccuracies  in  the  measurement  model.  By  using  a 
single  elemental  filter,  computational  loading  can  be  limited; 
additionally,  there  would  be  no  question  whether  the  multiple 
model  formulation  is  the  cause  of  this  behavior. 

Implementation  of  the  constant  turn  rate  target  acceleration 
model  (implemented  by  Suizu  [13]  in  the  extended  Kalman  filter 
formulation)  for  one  or  more  of  the  elemental  filters.  This 
model  more  closely  represents  high  performance  aircraft  evasive 
actions,  particularly  at  short  range. 

Change  hot  spot  separations  to  values  more  representative  of 
real  world  targets 

Add  pseudo-noise  to  the  position  and  velocity  states  to  account 
for  the  rotating  reference  frame  of  the  tracker 

Replace  the  current  impulsive  control  with  a  more  realistic 
control  response 


Improve  the  reacquisition  function  by  preventing  the  data  from 
the  diverging  filters  from  being  used  for  several  sample  periods 
following  reacquisition.  Currently,  if  the  MMAF  criteria 
(Bayesian  or  MAP)  selects  a  diverging  filter  following 
reacquisition,  it  is  included  even  though  its  estimate  is 
degraded.  By  preventing  the  use  of  these  estimates  for  a  few 
sample  periods,  the  conditional  probabilities  could  again 
converge  to  appropriate  values. 

Implement  the  algorithm  on  a  small  word-length  machine 

Investigate  the  implications  of  tracking  the  radiation  centroid, 
which  may  not  be  located  on  the  target,  rather  than  a  physical 
point  on  the  target.  Include  the  effects  of  the  hot-spot 
produced  by  laser  illumination  of  the  target  on  tracking 
performance. 
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Figure  A-2b.  Performance  Plot  for  TIFI-BASE 
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Figure  A-4a.  Performance  Plot  for  T2G 10MM-BASE 
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Figure  A-Sa.  Performance  Plot  for  T3G2M-BASE 
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Figure  A-9c.  Performance  Plot  for  T3O10MM-BASE 
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Figure  B-lf.  Performance  Plot  for  TIIIM-3AY 
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Figure  B-3d.  Performance  Plot  for  T2G2MM-BAY 
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i"ure  B-3f.  Performance  Plot  for  TZGZ'El-BAY 
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Performance  Plots  for  the  MAP  MMAF 
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Figure  C-lb.  Performance  Plot  for  TIMM-MAP 


Figure  C-2c.  Performance  Plot  for  T2G2MM-MAP 
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Figure  C-4a.  Performance 


Figure  C-4b.  Performance  Plot  for  T2G10MM-MAP  (Retuned) 
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Figure  C-Ac.  Performance  Plot  for  T2G10MM-MAP  (Retuned) 


Figure  C-4d.  Performance  Plot  for  T2G10MM-MAP  (Retuned) 
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Figure  C-7f.  Performance  Plot  for  T3G10MM 
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Figure  D-2b.  Performance  Plot  for  T2G2MM-XT 
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Figure  D-2c.  Performance  Plot  for  T2G2MM-XT 
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Figure  D-5f.  Performance  Plot  for  T3G10MM-XT 
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Figure  E-la.  Performance  Plot  for  TIMM-100HZ 
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Figure  E-ld.  Performance  Plot  for  TIMM-100HZ 
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Figure  E-2d.  Performance  Plot  for  T2G2MM-100HZ 
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Figure  E-3a.  Performance  Plot  for  T2G10MM-100HZ 
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Figure  E-4a.  Performance  Plot  for  T3G2MM-100HZ 
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As  noted  in  Section  6.2.2,  the  original  code  for  the  100  Hz  sample 
rate  resulted  in  loss  of  lock.  Subsequent  investigation  revealed  that 
the  reacquisition  revealed  that  the  reacquisition  subroutine  had  not 
been  implemented  correctly.  The  original  code  reset  the  covariance 
matrix  to  the  initial  values,  which  are  much  higher  than  steady  state 
values  and  have  no  cross-terms.  The  higher  than  appropriate  values 
caused  a  very  low  emphasis  on  incoming  data  at  a  point  when  the  measure¬ 
ments  should  be  weighted  heavily.  When  the  code  is  correctly  implemented, 
the  MMAF  does  not  lose  lock,  as  shown  in  the  following  plots  of  trajectory 
two  at  20  g's.  At  this  point,  there  is  no  explanation  for  the  significant 
bias  following  the  maneuver  recovery. 
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